AFRL-RY-WP-TR-2010-1024 


EFFECTS  OF  MULTIPLE  PHOTON  SCATTERING  IN 
DECIDUOUS  TREE  CANOPIES 

Michael  Greiner,  Bradley  D.  Duncan,  and  Matthew  P.  Dierking 

EO  Combat  ID  Branch 
EO  Sensor  Technology  Division 


DECEMBER  2009 
Interim  Report 


Approved  for  public  release;  distribution  unlimited. 

See  additional  restrictions  described  on  inside  pages 


STINFO  COPY 


AIR  FORCE  RESEARCH  LABORATORY 
SENSORS  DIRECTORATE 

WRIGHT-PATTERSON  AIR  FORCE  BASE,  OH  45433-7320 
AIR  FORCE  MATERIEL  COMMAND 
UNITED  STATES  AIR  FORCE 


NOTICE  AND  SIGNATURE  PAGE 


Using  Government  drawings,  specifications,  or  other  data  included  in  this  document  for  any 
purpose  other  than  Government  procurement  does  not  in  any  way  obligate  the  U.S.  Government. 
The  fact  that  the  Government  formulated  or  supplied  the  drawings,  specifications,  or  other  data 
does  not  license  the  holder  or  any  other  person  or  corporation;  or  convey  any  rights  or 
pennission  to  manufacture,  use,  or  sell  any  patented  invention  that  may  relate  to  them. 

This  report  was  cleared  for  public  release  by  the  USAF  88th  Air  Base  Wing  (88  ABW)  Public 
Affairs  Office  (PAO)  and  is  available  to  the  general  public,  including  foreign  nationals.  Copies  may 
be  obtained  from  the  Defense  Technical  Information  Center  (DTIC)  (http://www.dtic.mil). 

AFRL-RY -WP-TR-20 10-1 024  HAS  BEEN  REVIEWED  AND  IS  APPROVED  FOR 
PUBLICATION  IN  ACCORDANCE  WITH  ASSIGNED  DISTRIBUTION  STATEMENT. 


*//Signature// 


//Signature// 


LAWRENCE  J.  BARNES,  Project  Engineer  ROBERT  J.  FELDMANN,  Chief 


EO  Combat  ID  Branch 
EO  Sensor  Technology  Division 


EO  Combat  ID  Branch 
EO  Sensor  Technology  Division 


//Signature// 

BRIAN  C.  FORD,  Colonel,  USAF 
Chief,  EO  Sensor  Technology  Division 
Sensors  Directorate 


This  report  is  published  in  the  interest  of  scientific  and  technical  information  exchange,  and  its 
publication  does  not  constitute  the  Government’s  approval  or  disapproval  of  its  ideas  or  findings. 

*Disseminated  copies  will  show  “//Signature//”  stamped  or  typed  above  the  signature  blocks. 


REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including 
suggestions  for  reducing  this  burden,  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway,  Suite 

1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 

1.  REPORT  DATE  (DD-MM-YY) 

December  2009 

2.  REPORT  TYPE 

Interim 

3.  DATES  COVERED  (From  -  To) 

21  October  2007  -  20  December  2009 

4.  TITLE  AND  SUBTITLE 

EFFECTS  OF  MUFTIPFE  PHOTON  SCATTERING  IN  DECIDUOUS  TREE 
CANOPIES 

5a.  CONTRACT  NUMBER 

In-house 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

62204F 

6.  AUTHOR(S) 

Michael  Greiner,  Bradley  D.  Duncan,  and  Matthew  P.  Dierking 

5d.  PROJECT  NUMBER 

2003 

5e.  TASK  NUMBER 

11 

5f.  WORK  UNIT  NUMBER 

2003 112Y 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

EO  Combat  ID  Branch  (AFRL/RYJM) 

EO  Sensor  Technology  Division 

Air  Force  Research  Laboratory,  Sensors  Directorate 

Wright-Patterson  Air  Force  Base,  OH  45433-7320 

Air  Force  Materiel  Command,  United  States  Air  Force 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

AFRL-RY -WP-TR-20 1 0- 1 024 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory 

Sensors  Directorate 

Wright-Patterson  Air  Force  Base,  OH  45433-7320 

Air  Force  Materiel  Command 

United  States  Air  Force 

10.  SPONSORING/MONITORING 
AGENCY  ACRONYM(S) 

AFRL/RYJM 

11.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER(S) 

AFRL-RY- WP-TR-20 1 0- 1 024 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited. 

13.  SUPPLEMENTARY  NOTES 

PAO  Case  Number:  88ABW-09-4508;  Clearance  Date:  26  Oct  2009.  This  report  contains  color.  This  report  is  based  on  a 
PhD  dissertation  generated  to  fulfill  the  requirements  of  the  University  of  Dayton. 


14.  ABSTRACT 

Detecting  objects  hidden  beneath  forest  canopies  has  proven  to  be  a  difficult  task  for  optical  remote  sensing  systems. 
Rather  than  relying  upon  the  existence  of  gaps  between  the  leaves,  our  goal  was  to  use  the  light  that  is  scattered  from  the 
leaves  to  image  through  dense  foliage.  We  developed  a  Monte  Carlo  canopy  propagation  model  to  simulate  the  scattering 
of  light  through  a  maple  tree  canopy.  We  measured  several  forest  parameters,  including  the  gap  fraction  and  maximum 
leaf  area  density  of  a  real  test  canopy  and  applied  them  to  the  model.  We  ran  the  simulation  for  80°  illumination  and 
reported  on  the  results  in  the  ground  and  receiver  planes.  We  then  authenticated  the  validity  of  the  model  by  illuminating 
a  test  forest  at  an  80°  angle,  collecting  data  both  on  the  canopy  floor  and  in  a  monostatic  receiver,  and  comparing  the 
results  to  the  simulation.  Additionally,  we  examined  the  accuracy  of  the  model  in  accounting  for  seasonal  canopy 
variations  and  verify  the  simulation  with  experimental  results.  Lastly,  we  investigated  methods  for  boosting  the  signal-to- 
noise  ratio  (SNR)  of  detected  photons  and  make  SNR  calculations  for  various  illumination  angles. 


15.  SUBJECT  TERMS 

LADAR,  foliage  penetration,  foliage  scattering,  laser  propagation  modeling 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

OF  ABSTRACT: 

SAR 

18.  NUMBER 

OF  PAGES 

80 

19a.  NAME  OF  RESPONSIBLE  PERSON  (Monitor) 
Lawrence  J.  Barnes 

19b.  TELEPHONE  NUMBER  (Include  Area  Code) 

N/A 

a.  REPORT 

Unclassified 

b.  ABSTRACT 

Unclassified 

c.  THIS  PAGE 

Unclassified 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  Z39-18 


i 


Table  of  Contents 


Section  Page 


List  of  Figures . iv 

List  of  Tables . vi 

1.  Summary . 1 

2.  Introduction . 3 

2.1  Background . 3 

2.2  Leaf  Absorption  Spectrum . 3 

2.3  Methodology . 4 

3.  Methods,  Assumptions,  and  Procedures . 6 

3.1  Bidirectional  Scattering  Distribution  Functions . 6 

3.2  Monte  Carlo  Simulation . 9 

3.2.1.  Random  Propagation  Distance . 13 

3.2.2.  Leaf  Scattering  Angles . 17 

3.2.3.  Ground  and  Target  Scattering  Angles . 19 

3.2.4.  Test  Canopy  Parameters . 19 

3.3  Experimental  Verification . 21 

4.  Results  and  Discussions . 24 

4.1  Leaf  BSDF  Measurement  Results . 24 

4.1.1.  Calculation  of  Absorption  Coefficient . 24 

4.1.2.  Modeling  Specular  Reflection . 29 

4.1.3.  Modeling  Transmission  and  Diffuse  Reflection . 32 

4.1.4.  Leaf  Data  Interpolation . 33 

4.1.5.  Procedure  for  Constructing  the  BSDF . 35 

4.2  Monte  Carlo  Simulation  Results . 36 

4.3  Experimental  Verification  Results . 39 

4.4  Monte  Carlo  Results  at  Various  Illumination  Angles . 42 

4.4.1.  SNR  Calculations . 46 

4.4.2.  Results  of  Field  of  View  and  Range  Gate  Filtering . 48 

4.4.3.  Additional  SNR  Considerations . 50 

5.  Conclusions  and  Recommendations . 51 

6.  References . 53 

Appendix  A:  Calculation  of  Probability  of  a  Ballistic  Photon . 55 

Appendix  B:  APD  Field  of  View  Characterization . 59 

Appendix  C:  Test  Canopy  GPS  Waypoints . 63 

Appendix  D:  Cable  Delay  and  Dispersion . 65 

Appendix  E:  Selection  of  Random  Value . 67 

List  of  Acronyms,  Abbreviations  and  Symbols . 68 


iii 


List  of  Figures 


Figure  lAbsorption  Spectrum  of  Fresh  Poplar  Leaves13 . 4 

Figure  2  Photograph  of  the  BSDF  Measurement  Apparatus . 7 

Figure  3  Schematic  Diagram  of  the  BSDF  Measurement  apparatus . 7 

Figure  4  BSDF  of  (a)  a  Fresh  Common  Maple  Leaf,  and  (b)  an  Appreciably  Dried  Common  Maple  Leaf. . 9 

Figure  5  Example  Flemispheric  Image  Captured  Looking  Upward  within  a  Grove  of  Sugar  Maple  Trees . 9 

Figure  6  The  Canopy  is  Illuminated  by  a  Monostatic  Ladar  System  at  an  Angle  9inc . 10 

Figure  7  Monte  Carlo  Algorithm  Flow  Chart  Describing  Photon  Propagation  through  a  Tree  Canopy . 11 

Figure  8  Leaf  Orientation  Angles  (a)  are  Defined  in  the  (x,y,z)  Canopy  Coordinate  System,  while  Leaf  Scattering 

Angles  (b)  are  Defined  in  the  Leaf  Coordinate  System . 13 

Figure  9  Representative  Maple  Leaf  Area  Density  as  a  Function  of  Canopy  Height . 16 

Figure  10  Illustration  of  Region-to-Region  Propagation  in  a  Vertically  Segmented  Tree  Canopy . 17 

Figure  1 1  Example  BSDF  for  a  Sugar  Maple  Leaf  Illuminated  at  an  Incidence  Zenith  Angle  of  1 10° . 18 

Figure  12  Mean  Gap  Fraction  as  a  Function  of  Zenith  Angle  for  our  Nearby  Gove  of  Sugar  Maple  Trees . 21 

Figure  13  Image  of  the  Test  Canopy  as  Viewed  from  the  Tower . 22 

Figure  14  12x5  Measurement  Grid  Beneath  the  Canopy . 22 

Figure  15  Image  of  the  Path  Connecting  the  Entrance  Location  to  Target  as  Seen  from  the  Entrance  Location . 23 

Figure  16  Measured  BSDF  Data  of  Common  Maple  Leaves  for  Illumination  Angles  #7  of  (a)  0°  (b)  10°  (c)  20°  (d) 

30°  (e)  40°  (f)  50°  (g)  60°  (h)  70°  and  (i)  78° . 25 

Figure  17  Measured  BSDF  Data  of  Cottonwood  Leaves  for  Illumination  Angles  0 :  of  (a)  0°  (b)  10°  (c)  20°  (d)  30° 

(e)  40°  (f)  50°  (g)  60°  (h)  70°  and  (i)  78° . 26 

Figure  18  Comparison  of  the  Sugar  Maple  Leaf  BSDF  and  Spectralon®  BRDF  for  Normal  Illumination . 27 

Figure  19  Absorption  Coefficient  for  Maple  Leaves  as  a  Function  of  Illumination  Angle . 28 

Figure  20  Lambertian  Fit  to  the  BRDF  (left)  and  BTDF  (right)  for  (a)  Normal  Illumination  and  (b)  Illumination  at 

70° . 29 

Figure  21  Separation  of  the  Sugar  Maple  Leaf  Diffuse  and  Specular . 30 

Figure  22  Normalized,  Reverse  Rayleigh  Fits  of  the  Specular  Reflection  Data  for  Illumination  Angles  of  (a)  50°,  (b) 

60°,  (c)  70°,  and  (d)  78° . 31 

Figure  23  Fractional  Specular  Reflection  for  Maple  Leaves  as  a  Function  of  Incident  Angle . 32 

Figure  24  Fitted  Sugar  Maple  Leaf  BSDF  Curves  for  Illumination  at  70° . 32 

Figure  25  Sugar  Maple  Leaf  BRDF  Models  Created  by  Adding  the  Diffuse  and  Specular  Components  at  (a)  50°,  (b) 

60°,  (c)  70°,  and  (d)  78°  Illumination . 33 

Figure  26  Polynomial  Fits  of  the  (a)  0th,  (b)  1st,  and  (c)  2nd  Order  p-Coefficients  Describing  the  BTDF  Data  for 

Sugar  Maple  Leaves . 34 

Figure  27  Two  Dimensional  BRDF  Surface  Fit  for  Sugar  Maple  Leaves . 36 

Figure  28  Two  Dimensional  BTDF  Surface  Fit  for  Sugar  Maple  Leaves . 36 

Figure  29  Characteristic  Temporal  Waveform  from  those  Simulated  Photons  which  Strike  the  Canopy  Floor . 37 

Figure  30  Simulated  Ground  Plane  Beam  Footprint  for  an  Illumination  Angle  of  80° . 37 

Figure  31  Simulated  Beam  Footprint  Cross  Sections,  and  Best  Fit  Gaussian  Distributions,  in  the  Range  (a)  and 

Cross-Range  (b)  Dimensions . 38 

Figure  32  Spatial  Distribution  of  Simulated  RMS  Pulse  Widths  on  the  Virtual  Canopy  Floor  (in  ns) . 38 

Figure  33  Simulated  1-D  Temporal  PDF  Measured  in  Pupil  Plane  of  the  Virtual  Detector . 39 

Figure  34  Image  of  the  Avalanche  Photo-Diode  Setup . 40 

Figure  35  Characteristic  Temporal  Waveform  Measured  by  an  Upwards  Looking  Avalanche  Photodiode  Placed  on 

the  Canopy  Floor . 40 

Figure  36  Measured  Ground  Plane  Beam  Footprint  for  an  Illumination  Angle  of  80° . 41 

Figure  37  Measured  Beam  Footprint  Cross-Section  Data,  and  Best  Fit  Gaussian  Curves,  in  the  Range  (a)  and  Cross- 

Range  (b)  Dimensions . 41 

Figure  38  Spatial  Distribution  of  Actual  RMS  Pulse  Widths  Measured  on  the  Canopy  Floor . 42 

Figure  39  Actual  1-D  Temporal  PDF  Measured  with  a  Range  Gated  Intensified  CCD  Camera  Located  in  the  Tower 
. . 42 


IV 


Figure  40  Simulated  Beam  Footprint  on  the  Ground  for  0°  (a),  30  °  (b),  60  0  (c)  and  80  °  (d)  Canopy  Illumination 

Angles . 43 

Figure  41  Cross  Section  of  the  Simulated  Beam  Footprint  and  Best  Fit  Gaussian  Distribution  in  the  Range  (a)  and 

Cross-  Range  (b)  Dimensions . 43 

Figure  42  Photon  Arrival  Time  PDF’s  without  Range  Gating  or  Field  of  View  Filtering  for  10°  (a)  and  30°  (b) 

Canopy  Illumination . 45 

Figure  43  Scatter  Plot  of  Photon  Arrival  Time  versus  Angle  of  Arrival  onto  Detector  for  10°  (a)  and  30°  (b) . 46 

Figure  44  Contour  Plot  of  SNR  as  a  Function  of  Gate  Delay  and  Half  Field  of  View  after  Field  of  View  for  a  Gate 

Width  of  O.lps  for  10°  (a)  and  30°  (b)  Canopy  Illumination . 49 

Figure  45  Photon  Arrival  Time  PDF’s  after  both  Range  Gate  and  Field  of  View  Filtering  for  10°  (a)  and  30°  (b) 
Canopy  Illumination . 50 


v 


List  of  Tables 


Table  1  Sugar  Maple  Leaf  Reflection,  Transmission  and  Absorption  Coefficients  as  a  Function  of  Illumination 


Angle . 28 

Table  2  Eastern  Cottonwood  Leaf  Reflection,  Transmission  and  Absorption  Coefficients  as  a  Function  of 

Illumination  Angle . 28 

Table  3  Sugar  Maple  leaf  BTDF  q-Coefficients . 34 

Table  4  Sugar  Maple  leaf  BRDF  q-Coefficients . 34 

Table  5  Eastern  Cottonwood  leaf  BTDF  q-Coefficients . 35 

Table  6  Eastern  Cottonwood  leaf  BRDF  q-Coefficients . 35 

Table  7  Actual  Photon  Entrance  Locations  and  Simulated  Principle  Spot  Centers . 44 

Table  8  Best  Fit  Gaussian  Beam  Widths  and  RMS  Fitting  Errors . 44 

Table  9  Signal-to-Noise  Ratios,  Gate  Delays  and  Widths,  and  Half  Fields  of  View  for  Various  Canopy  Illumination 
Angles . 49 


Table  10  Signal-to-Noise  Ratios,  Gate  Delays  and  Widths,  and  Half  Fields  of  View  for  Various  Canopy  Illumination 


Angles 


50 


vi 


1.  Summary 


Detecting  objects  hidden  beneath  forest  canopies  has  proven  to  be  a  difficult  task  for  optical  remote 
sensing  systems.  Rather  than  relying  upon  the  existence  of  gaps  between  the  leaves,  our  goal  is  to  use  the 
light  that  is  scattered  from  the  leaves  to  image  through  dense  foliage.  We  develop  a  Monte  Carlo  canopy 
propagation  model  to  simulate  the  scattering  of  light  through  a  maple  tree  canopy.  We  measure  several 
forest  parameters,  including  the  gap  fraction  and  maximum  leaf  area  density  of  a  real  test  canopy  and 
apply  them  to  the  model.  We  run  the  simulation  for  80°  illumination  and  report  on  the  results  in  the 
ground  and  receiver  planes.  We  then  authenticate  the  validity  of  the  model  by  illuminating  a  test  forest  at 
an  80°  angle,  collecting  data  both  on  the  canopy  floor  and  in  a  monostatic  receiver,  and  comparing  the 
results  to  the  simulation.  Additionally,  we  examine  the  accuracy  of  the  model  in  accounting  for  seasonal 
canopy  variations  and  verify  the  simulation  with  experimental  results.  Lastly,  we  investigate  methods  for 
boosting  the  signal-to-noise  ratio  (SNR)  of  detected  photons  and  make  SNR  calculations  for  various 
illumination  angles. 

As  predicted  by  the  simulation,  we  experimentally  verified  that  a  large  number  of  photons  are  down 
scattered  just  after  entering  the  canopy,  creating  a  large  spot  on  the  ground  beneath  the  entrance  location. 
We  then  found  that  a  Gaussian  distribution  fit  the  range  and  cross-range  cross-sections  of  the  simulated 
and  measured  beam  footprints  with  great  accuracy.  We  also  found  that  the  standard  deviations  of  the  best- 
fit  Gaussian  distributions  were  on  the  same  order  in  both  simulation  and  experiment.  Additionally,  we 
observed  that  the  pulse  widths  of  the  waveforms  measured  on  the  ground  were  similar  in  both  shape  and 
magnitude  to  those  predicted  by  our  simulation.  Finally,  we  found  that  the  measured  1-D  temporal 
probability  density  function  (PDF)  of  photons  returning  to  the  receiver  closely  matched  the  PDF  predicted 
by  the  simulation. 

We  believe  that  any  mismatch  between  our  simulation  and  experimental  results  can  be  attributed  to  two 
factors.  First,  in  the  experiment  there  was  a  slight  misalignment  between  the  trajectory  of  the  initial  beam 
axis  and  the  path  cleared  in  the  tree  grove,  which  resulted  in  a  small  angular  error  in  the  pointing  of  our 
beam.  Second,  and  most  importantly,  the  simulation  results  represent  an  average  over  billions  of 
realizations  of  the  canopy  while  the  experiment  considers  only  one  real  forest  realization.  Considering 
these  factors,  we  believe  our  experimental  results  correlate  very  well  with  our  simulation,  leading  us  to 
conclude  that  our  model  is  valid. 

We  then  investigated  the  ability  of  the  canopy  propagation  model  to  handles  seasonal  canopy  variations. 
We  first  quantified  the  effects  of  seasonal  changes  on  two  forest  parameters:  the  bidirectional  scattering 
distribution  function  of  individual  leaves  and  the  maximum  leaf  area  density  of  the  canopy.  Then  by 
applying  these  parameters  to  our  Monte  Carlo  canopy  propagation  model  we  predicted  the  effect  of 
seasonal  changes  on  the  shape  and  size  of  the  beam  footprint  as  measured  on  the  ground.  We  found  that 
as  the  forest  health  declines  and  the  leaves  fall  from  the  trees  the  size  of  the  spot  on  the  ground  increases. 
We  then  experimentally  verified  the  results  of  the  model  by  illuminating  a  real  forest  with  a  780nm  beam 
and  sampling  the  beam  footprint  with  an  upwards  facing  avalanche  photo-diode  (APD)  beneath  the 
canopy.  The  same  trend  was  seen  in  the  experimentally  measured  data  as  was  predicted  by  the  simulation. 

Finally,  we  used  our  Monte  Carlo  canopy  propagation  model  to  make  predictions  about  detecting  objects 
through  the  canopy  at  several  different  illumination  angles.  We  have  used  our  Monte  Carlo  propagation 
model  to  simulate  the  scattering  of  photons  through  a  dense  canopy  using  several  different  illumination 
angles.  We  found  that  the  scattered  beam  footprint  on  the  canopy  floor  is  primarily  located  beneath  the 
entrance  location  of  photons  into  the  canopy  for  all  illumination  angles.  As  a  result,  illuminating  the 
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canopy  above  an  anticipated  target,  rather  than  aiming  directly  at  it,  places  more  photons  on  the  target  and 
immensely  increases  the  signal-to-noise  ratio  of  detectable  returning  photons. 

Additionally,  we  examined  the  effects  of  using  a  narrow  field  of  view,  range  gated  camera.  We  found  that 
multiply-scattered  photons  spend  more  time  in  the  canopy  than  do  singly-scattered  photons,  and  therefore 
generally  arrive  at  the  virtual  receiver  plane  at  a  later  time.  The  use  of  a  range  gated  camera  allows  us  to 
look  past  the  initial  return  of  backscattered  noise  photons,  increasing  the  signal-to-noise  ratio.  Multiple 
canopy  scattering  also  causes  spatial  and  angular  dispersion  of  photons  throughout  the  canopy.  Thus, 
photons  arriving  at  later  times  most  likely  also  arrive  at  larger  angles  which  allows  many  noise  photons  to 
be  filtered  using  a  narrow  field  of  view  camera. 

We  found  that  by  illuminating  directly  above  the  desired  target  and  using  a  narrow  field  of  view,  range 
gated  camera  we  found  that  we  can  boost  the  signal-to-noise  ratio  on  the  order  of  15dB.  While  our  results 
are  specific  to  the  exact  geometry  and  canopy  architecture  used  in  this  simulation,  we  have  demonstrated 
in  principle  that  we  may  boost  the  signal-to-noise  ratio  of  detected  multiply  scattered  photons  using  a 
narrow  field  of  view,  range  gated  camera,  ultimately  providing  the  possibility  of  imaging  obscured  targets 
embedded  within  a  dense  forest  canopy  at  low  illumination  angles. 


2 


2.  Introduction 


2.1  Background 

An  on-going  problem  for  the  military  is  finding  and  identifying  objects  concealed  by  foliage  or 
camouflage.  Forest  canopies  are  a  natural  barrier  to  air  or  satellite  detection  using  both  optical  and 
traditional  radar  techniques.  The  traditional  approach  to  imaging  objects  hidden  by  foliage  relies  on 
synthetic  aperture  radar  (SAR)  to  penetrate  through  the  leaves. 1,2,3  Because  the  wavelength  is  typically  on 
the  order  of  several  meters,  the  light  passes  through  the  vegetation  with  minimal  attenuation  and 
therefore,  this  technique  has  proven  to  be  successful  at  detecting  targets.  However,  also  because  the 
wavelength  is  on  the  order  of  several  meters,  the  resolution  is  limited  and  a  clear  high  definition  image 
cannot  be  formed.  Thus,  in  the  interest  of  achieving  high-fidelity  target  identification  and  minimizing 
false  alarms,  researchers  have  been  exploring  the  advantages  of  active  laser  imaging  since  the  early 
1980s.4,5 

Current  optical  foliage  penetration  methods  depend  heavily  upon  looking  through  gaps  between  the 
leaves  and  branches.1  For  example,  as  an  airborne  sensor  is  flown  over  a  forested  area,  small  pieces  of  a 
target  under  the  canopy  may  be  sequentially  imaged  and  later  pieced  together  to  form  a  composite 
image.6,7  Unfortunately,  the  gap  fraction  quickly  drops  to  zero  as  the  foliage  density  and/or  zenith  angle 
increases,  leaving  very  few  opportunities  to  exploit  gaps  in  the  canopy.8  Rather  than  relying  upon  the 
existence  of  gaps  between  the  leaves,  our  goal  is  to  use  the  light  that  is  scattered  from  the  leaves  to  image 
through  dense  foliage. 

2.2  Leaf  Absorption  Spectrum 

At  optical  wavelengths,  absorption  occurs  within  leaves  at  specific  wavelengths  due  to  chlorophyll,  water, 
and  other  leaf  materials.9  Light  absorbed  by  chlorophyll  may  be  converted  into  energy  through 
photosynthesis.10  Chlorophyll  has  high  absoiption  in  the  visible  spectrum,  ranging  from  the  blue  (around 
445  nm)  to  the  red,  around  645nm.  Leaf  water  is  also  a  major  constituent  of  fresh  leaves,  representing 
around  66%  of  a  leafs  weight.11  Leaf  water  absoiption  has  large  absoiption  features  in  the  infrared  range, 
beginning  around  1400nm,  but  close  to  zero  absoiptance  in  the  visible  or  near  infrared  part  of  the 
spectrum.12  The  overall  leaf  absoiption  spectrum  is  shaped  by  the  spectra  of  chlorophyll  and  water,  as 
well  as  its  many  other  components.  An  example  of  the  absoiption  spectrum  of  a  fresh  poplar  leaf  can  be 
seen  in  Figure  l.13 
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Figure  1  Absorption  Spectrum  of  Fresh  Poplar  Leaves 


Light  striking  any  object  must  be  reflected,  transmitted,  or  absorbed.  Thus,  the  figure  is  broken  into  three 
distinct  regions:  transmittance  T ,  reflectance  R  ,  and  absorptance  A  ,  where  the  absoiptance  is 
determined  from  T  and  R  through  the  relationship  A  =  1  —  R  —  T  .  For  the  majority  of  foliage  there  is  a 
band  of  the  spectrum,  ranging  from  800nm  to  1400nm,  where  the  absoiption  of  light  dips  sharply,  and 
transmission  and  reflection  are  relatively  large  and  nearly  equal.  A  wavelength  found  in  this  low 
absorption  band  could  potentially  be  used  to  penetrate  the  leaves  and,  in  effect,  see  through  them  to 
reflective  objects  behind  them.  These  curves  are  typical  of  most  types  of  deciduous  tree  leaves.14 


2.3  Methodology 

In  Section  2  we  discuss  several  key  concepts  related  to  this  project.  We  examine  the  spectral  properties  of 
individual  deciduous  tree  leaves,  as  well  as  the  bidirectional  scattering  distribution  functions  of  sugar 
maple  and  eastern  cottonwood  leaves  for  1064nm  light.  We  also  discuss  the  forward  and  reverse 
transformations  of  random  variables,  which  we  use  to  select  random  distances  and  angles  in  our  Monte 
Carlo  algorithm.  Finally,  we  take  a  look  at  the  fundamental  principles  of  range  gated  imaging  and  define 
several  key  terms. 

In  Section  3  we  introduce  a  canopy  propagation  model  which  will  allow  us  to  simulate  the  propagation 
and  scattering  of  photons  through  a  forest  canopy  with  one  hundred  percent  foliage  obscuration.  We 
describe,  in  detail,  the  Monte  Carlo  algorithm  along  with  the  individual  probability  density  functions 
governing  each  of  the  random  variables  we  use  in  this  simulation.  We  discuss  the  measurement  of  several 
parameters  from  a  nearby  tree  grove  and  apply  them  to  our  model.  By  using  real  forest  parameters,  we 
intend  to  achieve  the  most  accurate  and  reproducible  results.  We  report  the  findings  of  our  simulations 
and  investigate  the  validity  of  the  model  by  comparing  several  predicted  phenomena  to  that  which  we 
measured  by  propagating  a  real  beam  through  the  canopy  of  a  local  maple  tree  grove. 

In  Section  4  we  investigate  the  application  of  the  model  to  seasonal  changes  in  the  canopy.  As  the  water 
in  the  leaves  is  replaced  with  air  the  index  of  refraction  lowers  and  the  scattering  properties  change.  We 
therefore  look  at  the  seasonal  dependence  of  the  BSDF  for  individual  sugar  maple  and  eastern 
cottonwood  leaves.  We  also  examine  the  seasonal  dependence  of  the  leaf  number  density  within  a  foliated 
forest  by  capturing  hemispheric  images  at  a  set  grid  of  locations  beneath  the  canopy.  We  apply  the 
measured  canopy  parameters  to  our  Monte  Carlo  model  and  report  on  the  simulation  results.  We  then 
illuminate  our  test  canopy  with  a  real  beam  and  perform  waveform  measurement  experiments  to  validate 
the  accuracy  of  the  model’s  prediction  to  seasonal  changes. 
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In  Section  5  we  examine  the  possibility  of  using  scattered  light  for  imaging  through  a  healthy,  foliated 
canopy.  Under  the  assumption  that  the  canopy  propagation  model  is  valid,  we  run  the  simulation  for  a 
variety  of  illumination  angles  and  make  predictions  on  the  prospect  of  imaging.  We  derive  an  expression 
for  the  signal  to  noise  ratio  of  the  detection  process  as  a  function  of  the  number  of  signal  and  noise 
photons  collected  by  the  detector.  We  then  examine  the  effects  of  using  a  narrow  field  of  view  camera  as 
well  as  implementing  range  gated  imaging  and  report  on  their  expected  boost  in  signal  to  noise  ratio. 
Finally,  we  present  a  summary  of  our  work  and  our  key  conclusions  in  Section  6. 
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3.  Methods,  Assumptions,  and  Procedures 


3.1  Bidirectional  Scattering  Distribution  Functions 

The  scattering  and  transmission  properties  of  trees  and  forests  are  topics  of  great  interest  in  the  field  of 
optical  remote  sensing.  As  leaves  dominate  the  scattering  of  light  within  a  forest  canopy,  knowledge  of 
their  optical  properties  is  essential  in  order  to  model  the  interaction  of  photons  with  the  canopy.  Optical 
properties  of  leaves  have  been  the  subject  of  many  studies  in  recent  years.  These  studies  have  found 
applications  in  several  fields,  including  photobiology,  agriculture  and  remote  sensing. 12,15,16 

This  section  focuses  upon  the  bidirectional  scattering  distribution  function  (BSDF)  of  maple  and 
cottonwood  leaves.  The  BSDF  of  a  surface  is  the  ratio  of  the  scattered  radiance  to  incident  irradiance  at  a 
given  wavelength.  The  function  is  dependent  upon  two  directional  angles,  the  angle  of  illumination  and 
the  angle  at  which  light  is  scattered.  The  BSDF  is  typically  split  into  reflected  and  transmitted 
components,  which  are  treated  separately  as  the  bidirectional  reflectance  distribution  function  (BRDF) 
and  the  bidirectional  transmittance  distribution  function  (BTDF),  respectively. 

Few  studies  have  been  performed  on  the  scattering  functions  of  individual  leaves  because  it  has  typically 
been  assumed  that  these  functions  are  simply  Lambertian  in  nature.17,18  However,  measurements  have  not 
always  supported  this  assumption.  In  fact,  some  studies,  including  our  own,  have  found  a  strong  specular 
reflection  component  accompanying  the  diffuse  Lambertian  scattering  for  higher  illumination  angles.19,20 

It  has  been  demonstrated  that  the  wavelength  of  the  illuminating  laser  dramatically  changes  the  scattering 
properties  of  leaves.  Figure  1  for  example,  depicts  the  reflectance  and  transmittance  spectra  of  fresh 
poplar  leaves.  Although  poplar  ( Populus  canadensis )  is  not  one  of  the  species  investigated  in  this 
research,  the  general  trend  in  its  absoiption  spectrum  is  typical  of  all  deciduous  leaves  and  is  therefore 
relevant  to  our  work.  Figure  1  is  broken  into  three  distinct  regions:  transmittance  T ,  reflectance  R  ,  and 
absoiptance  A  ,  where  the  absoiptance  is  determined  from  T  and  R  through  the  relationship 
A  =  1  —  R  —  T  .  Notice  that  the  visible  region  is  characterized  by  high  absoiptance.  There  are  also  strong 
absoiption  peaks  in  the  infrared.  However,  there  is  a  region  in  the  near-IR,  between  800  and  1300  nm, 
where  absorptance  is  minimized.  Selecting  a  wavelength  in  this  region  therefore  provides  the  largest 
amount  of  light  reflected  by  and  transmitted  through  the  leaves.  It  will  thus  be  the  goal  of  this  research  to 
characterize  the  BSDFs  of  two  deciduous  leaf  species  in  the  local  Dayton,  OH  area  (i.e.,  Sugar  Maple  and 
Eastern  Cottonwood)  at  a  single,  near-IR  wavelength. Measurement  of  BSDF  data  from  individual  leaves 
is  performed  through  the  use  of  the  goniometric  apparatus  shown  in  Figure  2  and  depicted  schematically 
in  Figure  3.  A  linearly  polarized  1064  nanometer  pulsed  laser  (>6pJ  pulse  energy)  is  directed  along  the 
axis  of  a  stationary  optical  rail.  A  beam  splitter  is  used  to  send  half  of  the  beam  to  the  leaf  for  scattering 
and  the  other  half  to  an  energy  level  detector.  The  beam  directed  towards  the  energy  detector  is  first 
incident  upon  a  50%  reflective  Spectralon®  disc.  The  measured  radiant  energy  reflected  from  this  disc  is 
simply  used  to  monitor  pulse-to-pulse  energy  fluctuations  so  that  data  can  later  be  normalized  with 
respect  to  variations  in  pulse  energy. 

Two  neutral  density  (ND)  filters,  mounted  on  the  rail  after  the  beam  splitter,  are  used  to  attenuate  the 
power  of  the  laser  in  order  to  avoid  damaging  the  transmission  and  reflection  detectors.  The  ND  filters  are 
also  tilted  slightly  in  order  to  avoid  direct  reflections  into  the  beam  path.  Notice,  though,  that  any 
deflection  of  the  beam  due  to  the  first  filter  is  counteracted  by  the  second  so  that  the  path  of  the  beam 
remains  along  the  rail  axis. 
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Figure  2  Photograph  of  the  BSDF  Measurement  Apparatus 


After  attenuation,  the  laser  beam  is  directed  onto  the  leaf,  which  is  mounted  in  a  goniometer  for 
measurement.  The  goniometer  has  two  separate,  coaxial  rotation  stages  that  are  motor  driven  and 
independently  controlled  by  a  computer  which  also  tracks  the  leafs  rotational  position  relative  to  the 
beam  path.  A  second  optical  rail  is  mounted  on  one  of  the  rotation  stages  with  reflection  and  transmission 
detectors  equidistant  from  the  leaf  at  opposite  ends  of  the  rail.  The  active  region  of  the  detectors  is  small 
enough  to  give  point  measurements  and  avoid  angular  averaging  of  the  data.  The  leaf  is  mounted  above 
this  rail  on  a  second  rotation  stage  which  is  driven  by  a  separate  motor.  In  this  way  the  illumination  angle 
of  the  beam  on  the  leaf  can  be  adjusted  independently  of  the  observation  angles. 

Each  rotation  stage  is  capable  of  360°  of  rotation,  allowing  the  reflection  and  transmission  to  be  observed 
at  any  angle.  The  rail  rotation  angle  0D  is  considered  to  be  at  0°  when  the  transmission  detector  is  at  its 
leftmost  position  and  would  see  the  beam  passing  through  the  leaf  with  no  deflection.  The  leaf  rotation 
angle  0 1  is  considered  to  be  at  0°  when  the  leaf  surface  is  normal  to  the  incident  beam. 


Transmission 


Figure  3  Schematic  Diagram  of  the  BSDF  Measurement  apparatus 
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The  motors  of  both  the  detector  rail  and  leaf  mount  are  controlled  by  an  automated  LabView®  program 
which  rotates  the  detectors  and  leaf  to  predetermined  angles.  At  each  set  of  angles  the  program  takes  256 
samples  from  each  detector.  The  samples  are  then  averaged,  and  a  mean  value  is  stored  for  each  detector. 
The  program  then  moves  the  leaf  and/or  detectors  to  the  next  set  of  angles,  and  the  process  is  repeated 
until  all  the  desired  angle  permutations  are  exhausted.  Additionally,  the  entire  procedure  is  repeated  for 
four  different  leaves.  The  four  mean  values  for  each  detector  are  then  averaged  into  a  single  value  for 
every  angle  permutation.  The  result  is  the  average  of  1024  samples  taken  from  four  different  leaves.  The 
additional  data  samples  are  collected  in  order  to  reduce  any  leaf-to-leaf  variance  in  the  final  readings. 

The  diameter  of  the  Gaussian  beam  is  optically  set  to  approximately  1  cm  to  provide  a  large  illumination 
footprint.  This  is  done  for  two  reasons.  First,  laser  speckle  is  reduced  by  increasing  the  size  of  the 
illumination  beam.  Reducing  laser  speckle  in  turn  reduces  variance  in  readings  made  by  the  detectors. 
Second,  because  a  leafs  structure  is  not  homogeneous  (i.e.,  leaves  have  veins  and  stems,  as  well  as  other 
fine  cellular  structures  running  through  them),  it  is  important  to  illuminate  a  large  enough  area  of  the  leaf 
to  ensure  good  spatial  averaging. 

Speckle  reduction  and  spatial  averaging  are  also  enhanced  by  periodically  translating  the  leaf  within  the 
path  of  the  beam.  In  addition  to  being  mounted  on  a  rotation  stage,  the  leaf  is  mounted  on  a  vertical 
translation  stage,  shown  in  Figure  2,  which  is  driven  at  a  frequency  of  approximately  1  Hz.  Translating 
the  leaf  up  and  down,  coupled  with  the  use  of  a  broad  illumination  beam,  creates  a  large  effective  area 
over  which  the  leaf  is  sampled.  As  a  result,  localized  leaf  structure  has  minimal  effect  on  the 
measurements. 

BSDF  data  was  collected  from  Sugar  Maple  and  Eastern  Cottonwood  leaves  found  in  the  Dayton,  Ohio 
area  during  the  weeks  of  May  21,  2006  and  May  28,  2006.  For  each  leaf  species  both  the  bidirectional 
reflectance  (BRDF)  and  transmittance  (BTDF)  distribution  functions  were  measured  for  incident  angles 
of  0,  10,  20,  30,  40,  50,  60,  70,  and  78°.  Illumination  angles  higher  than  78°  were  not  used  because  the 
projected  beam  waist  becomes  larger  than  the  physical  width  of  the  leaf  as  it  rotates  to  higher  angles.  The 
scan  of  a  single  illumination  angle  then  involves  taking  measurements  of  the  scattered  radiance  at  many 
detector  angles  about  the  leaf  surface.  We  collected  data  at  detector  angles  spanning  a  range  from  -85°  to 
+85°  about  the  leaf  surface  normal,  in  increments  of  5°. 

Because  scattering  properties  change  as  leaves  dry  out,  scan  durations  must  be  kept  short.9  Preliminary 
experiments  showed  a  strong  correlation  between  the  freshness  of  a  leaf  and  its  scattering  characteristics. 
For  example,  Figure  4  shows  the  BSDF  of  a  (a)  fresh  Sugar  Maple  leaf  illuminated  with  normally 
incident  light,  and  (b)  the  same  leaf  left  out  to  dry  overnight.  In  these  figures,  light  incident  from  the  left 
(illustrated  by  the  arrow)  illuminates  a  leaf  (bold  line)  at  normal  incidence.  The  resulting  transmission 
profile  (right  hand  lobes)  and  reflection  profile  (left  hand  lobes)  are  plotted  in  polar  coordinates.  There  are 
appreciable  effects  due  to  leaf  drying,  as  can  be  seen  in  the  difference  in  the  size  of  the  reflection  and 
transmission  lobes  between  the  two  plots.  In  this  case,  as  water  in  the  leaf  is  replaced  with  air  we 
observed  that  the  reflectance  of  leaves  increases  while  the  transmittance  decreases.  Notice  that  the  shape 
of  the  two  lobes  remains  relatively  constant,  though  the  area  encompassed  by  each  changes  appreciably. 

Because  of  the  effect  of  leaf  drying  on  the  scattering  profile,  it  is  necessary  to  ensure  that  the  leaves 
remain  fresh  during  the  entire  measurement  procedure.  In  order  to  combat  the  effects  of  leaf  drying,  a 
single  leaf  was  used  to  scan  only  three  illumination  angles  before  it  was  replaced  with  a  new,  fresh  leaf. 

As  the  scan  of  a  single  illumination  angle  takes  approximately  25  minutes,  limiting  the  use  of  a  single  leaf 
to  less  than  90  minutes  proved  short  enough  that  drying  effects  were  not  noticeable  in  any  leaves,  or  our 
final  data. 


12 


0 


Figure  4  BSDF  of  (a)  a  Fresh  Common  Maple  Leaf,  and  (b)  an  Appreciably  Dried  Common  Maple  Leaf. 


3.2  Monte  Carlo  Simulation 

Detecting  objects  beneath  forest  canopies  is  a  difficult  task  for  optical  remote  sensing  systems.  Current 
foliage  penetration  methods  depend  heavily  upon  the  ability  to  look  through  gaps  between  the  leaves  and 
branches. 1  For  example,  as  an  airborne  sensor  is  flown  over  a  forested  area,  small  pieces  of  a  target  under 
the  canopy  may  be  sequentially  imaged  and  later  pieced  together  to  form  a  composite  image.21’22 
Unfortunately,  as  shown  in  Figure  24,  the  gap  fraction  quickly  drops  to  zero  as  the  foliage  density  and/or 
zenith  angle  increases,  leaving  very  few  opportunities  to  exploit  gaps  in  the  canopy.22  Developing  the 
capability  to  look  directly  through  leaves  while  not  relying  on  the  existence  of  gaps  will  greatly  aid  in  the 
imaging  of  objects  beneath  the  canopy. 

Imaging  from  scattered  photons  presents  new  and  unique  challenges.  The  scattering  distributions  of 
individual  leaves  have  been  studied  extensively,  but  understanding  what  happens  to  light  as  it  is  scattered 
through  a  random  collection  of  leaves  is  a  more  complex  problem.12’24  In  particular,  after  propagating 
through  a  canopy  a  light  beam  will  experience  spatial,  angular,  and  temporal  dispersion.  The  distribution 
of  photons  exiting  the  canopy  will  then  be  physically  wider,  de-collimated,  and  temporally  dispersed  from 
the  incident  beam.25 


Figure  5  Example  Flemispheric  Image  Captured  Looking  Upward  within  a  Grove  of  Sugar  Maple  Trees 
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The  propagation  of  light  through  a  forest  canopy  was  modeled  by  means  of  an  extensive  Monte  Carlo 
simulation  in  which  multiple  scattering  events  within  the  canopy  are  treated  as  a  sequence  of  interactions 
between  a  single  photon  and  a  discrete  scatterer.  Simplification  lies  in  the  fact  that  the  end  result  does  not 
come  through  looking  at  the  canopy  as  a  whole,  but  rather  by  considering  each  photon  individually  and 
each  interaction  sequentially.  The  advantage  of  using  Monte  Carlo  methods  is  that  multiple  scattering 
events  can  be  calculated  without  complex  analysis,  as  only  single  variable  scattering  probability  density 
functions  are  required. 

We  use  our  Monte  Carlo  algorithm  to  track  the  propagation  of  photons  through  a  simulated  canopy  which 
we  have  modeled  as  an  elliptically  cylindrical  collection  of  leaves  that  are  randomly,  but  not  uniformly, 
distributed.  As  shown  in  Figure  25,  photons  are  launched  at  an  illumination  angle  of  0iNC  from  an 
elevated  monostatic  direct  detection  ladar  system  and  are  aimed  at  an  assumed  primary  target  located  on 
the  ground  at  the  geometric  center  of  the  canopy.  We  have  also  assumed  a  secondary  ground  target 
beneath  the  location  where  photons  first  enter  the  canopy,  as  preliminary  simulations  showed  a  large  spot 
on  the  ground  at  this  location  due  to  initial  down-scattering.  Both  virtual  targets  are  1  Om  in  diameter  and 
are  tilted  so  that  their  surface  normals  are  pointed  back  toward  the  receiver. 

The  algorithm  used  is  based  on  the  cloud  propagation  model  presented  by  E.  Bucher  and  is  depicted  in  the 
flow  chart  provided  in  Figure  26.25  According  to  our  model,  a  single  incident  photon  enters  the  canopy 
and  travels  a  random  distance  before  it  experiences  its  first  scattering  interaction.  On  this  initial 
propagation  the  photon  cannot  stray  from  its  trajectory  and  therefore  may  only  encounter  a  leaf  or  the 
primary  target  at  which  it  is  initially  aimed.  However,  there  are  typically  four  possible  options  after 
random  photon  propagation:  (i)  leaf  interaction,  (ii)  target  interaction,  (iii)  ground  interaction,  or  (iv) 
canopy  departure. 

If  the  photon  remains  within  the  bounds  of  the  canopy  after  its  initial  propagation  event,  then  it  must  have 
struck  a  leaf,  and  either  absoiption  or  scattering  will  occur.  We  determine  whether  or  not  the  photon  is 
absorbed  by  selecting  a  value  of  a  random  variable  uniformly  distributed  between  0  and  1 .  If  the  random 
value  is  smaller  than  the  leaf  absoiption  coefficient,  then  the  photon  is  absorbed  and  we  launch  a  new 
photon  into  the  canopy.  Otherwise,  the  photon  is  scattered  from  the  leaf,  in  which  case  we  select  random 
leaf  orientation  angles  (0 \,  and  <j)\)  and  random  leaf  scattering  angles  (0S  and  <j)$). 


Figure  6  The  Canopy  is  Illuminated  by  a  Monostatic  Ladar  System  at  an  Angle  0INc 
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Figure  7  Monte  Carlo  Algorithm  Flow  Chart  Describing  Photon  Propagation  through  a  Tree  Canopy 


As  shown  in  Figure  27  (a),  the  leaf  orientation  angles  ( 0  L  and  (f>\)  are  defined  in  the  canopy  coordinate 
system,  where  is  the  zenith  angle  made  between  the  z-axis  and  the  leaf  normal  vector  n  ,  and  <j)\  is  the 
azimuth  angle  between  the  x-axis  and  the  projection  of  the  leaf  normal  vector  onto  the  x-y  plane.  The 
canopy  coordinate  system  is  in  turn  defined  so  that  the  origin  is  located  on  the  ground  at  the  center  of  the 
primary  target,  where  the  positive  x-axis  is  parallel  to  the  ground  and  points  directly  away  from  the 
transmitter,  the  z-axis  points  straight  up,  and  the  y-axis  is  defined  according  to  the  right-hand  rule.  As 
shown  in  Figure  27  (b),  the  leaf  scattering  angles  (Os  and  <f>s)  are  then  defined  in  the  leaf  coordinate 
system,  where  Os  is  the  zenith  angle  made  between  the  leaf  surface  normal  and  the  scattered  photon’s 
unit  propagation  vector  s  ,  and  (j>s  is  the  azimuth  angle  made  between  the  projections  of  s  and  the 
photon’s  unit  propagation  vector  p  as  it  is  incident  upon  the  leaf  surface. 

The  photon  propagation  angles  (0 and  (j>),  defined  in  the  canopy  coordinate  system  as  shown  in  Figure  27 
(a),  are  then  determined  using  a  simple  matrix  transformation.  In  particular,  after  scattering,  the  scalar 
components  of  the  unit  vector  describing  the  propagation  direction  s  can  be  expressed  in  terms  of  the 
leaf  coordinate  system  as 


sin^  cos  (j)s 

ys 

= 

sin  0S  sin  </>s 

zs\ 

cos  ds 
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This  vector  must  then  be  transformed  into  the  canopy  coordinate  system  by  applying  rotation  angles 
and  9  L  about  the  z-  and  y-axes,  respectively,  according  to 


X 

cos  6,  cos  <j>L 

~  sin 

sin  9l  cos  (j)L 

sin  9S  cos  <f>s 

y 

= 

cos  6l  sin  (j)L 

cos^L 

sin  9L  sin  (f>L 

sin  9S  sin  (fts 

z 

-  sin  9, 

0 

cos  9t 

cos  9S 

From  Eq.  (3.2)  the  desired  photon  propagation  angles  can  then  be  determined: 


9  =  tan 


-l 


V 


2  ,  2 

x  +  y 


and 


(f)  =  tan 


(A 

UJ 


(3) 


Next,  if  the  photon  reaches  the  ground  plane,  its  location  is  compared  to  the  locations  of  the  two  targets.  If 
the  photon  strikes  either  of  the  targets,  both  of  which  are  assumed  to  have  unit  reflectance,  we  save  the 
spatial,  temporal  and  angular  ground  plane  data  and  select  random  target  scattering  angles  ( 0T  and  <j)j). 
The  target  scattering  angles  are  then  defined  in  a  target  coordinate  system  analogous  to  the  leaf  coordinate 
system,  in  the  sense  that  the  zenith  and  azimuth  angles  are  determined  by  the  normal  vector  of  the  target. 
Also,  by  substituting  ( 0  y  and  c/)\)  for  ( 0  s  and  <j>s),  the  same  matrix  transformation  described  by  Eqns.  (1)- 
(3)  can  be  applied  to  determine  the  new  photon  propagation  angles  in  the  canopy  coordinate  system. 

If,  however,  the  photon  reaches  the  ground  plane  but  does  not  interact  with  a  target,  it  must  have  struck 
the  ground,  which  we  also  assume  to  have  unit  reflectance.  In  this  case  we  save  the  spatial,  temporal,  and 
angular  ground  plane  data  and  select  random  ground  scattering  angles  (9 g  and  </>(,).  Note  that  the  ground 
scattering  angles  are  defined  in  the  canopy  coordinate  system  and  therefore  describe  the  new  photon 
propagation  angles  without  further  transformation. 

Finally,  if  the  random  propagation  distance  places  the  photon  outside  the  bounds  of  the  canopy  we 
propagate  the  photon  back  to  the  receiver  plane.  To  match  our  experimental  conditions,  the  receiver  plane 
in  our  model  was  assumed  to  be  located  45  m  above  the  ground  at  a  range  of  255  m  (for  80°  illumination) 
and  tilted  such  that  its  surface  normal  is  directed  toward  the  primary  target  at  the  center  of  the  canopy.  If  a 
photon  hits  the  receiver  plane  within  a  25  cm  diameter  pupil  centered  on  the  monostatic  transmit/receive 
ladar  axis  and  within  a  45°  field  of  view,  then  we  save  the  spatial,  temporal,  and  angular  receiver  plane 
data  and  launch  a  new  photon.  Otherwise  the  photon  is  considered  undetected,  we  do  not  save  any  data 
and  then  we  launch  a  new  photon.  In  our  model,  photons  that  exit  the  canopy  traveling  nominally  away 
from  the  receiver  plane  will  arrive  at  the  detector  in  the  negative  time  domain.  These  photons  are  also 
discarded  as  undetected. 


12 


z 


Figure  8  Leaf  Orientation  Angles  (a)  are  Defined  in  the  (x,y,z)  Canopy  Coordinate  System,  while  Leaf 
Scattering  Angles  (b)  are  Defined  in  the  Leaf  Coordinate  System 

Each  photon  can  travel  any  distance  in  any  direction  according  to  a  number  of  random  variables  which 
must  be  examined  at  every  step.  In  particular,  the  canopy  is  described  completely  by  the  probability 
density  functions  describing  the  following  random  variables:  the  propagation  distance  d  ;  the  probability 
that  the  photon  is  reflected  R,  transmitted  T,  or  absorbed  A  by  the  leaf;  the  leaf  angular  orientation  angles 
0  l  and  the  leaf  scattering  angles  0  s  and  </>s',  the  ground  scattering  angles  0(,  and  </>(,',  and  the  target 
scattering  angles  0j  and  (f>\.  The  probability  density  functions  for  each  of  the  random  variables  will  be 
described  in  the  following  sub-sections.  For  convenience,  a  simple  method  for  generating  specific  values 
of  a  random  variable,  given  its  probability  density  function,  is  discussed  in  the  appendices. 

3.2.1.  Random  Propagation  Distance 

For  a  homogeneous  medium  of  constant  leaf  number  density,  the  distance  d  that  the  photon  travels 
between  scattering  events  is  a  random  variable  described  by  the  following  exponential  probability  density 
function, 25 


Pj(d) 


1  d 

—  exp - 

D  D 


(4) 


where  the  parameter  D  is  the  mean  free  path,  or  average  distance  traveled  by  the  photon  between 
interactions.  The  mean  free  path  D(z,  0)  is  in  turn  inversely  proportional  to  the  product  of  the  mean 

projected  leaf  area  Ap  (0)  in  the  direction  of  photon  propagation  and  the  leaf  number  density  ,/V(z)  of 
the  canopy  in  the  region  surrounding  the  photon.  In  particular, 


D(z,0) 


1 

ap(9\n(z) 


(5) 


where  0  is  defined  in  Eq.  (3)  and  z  is  the  vertical  height  of  the  photon  within  the  canopy.26 
3.2.1  A.  Mean  Projected  Leaf  Area 

Because  leaves  mostly  face  upwards  towards  the  sun,  photons  traveling  at  zenith  angles  f?far  off  the 
vertical  will  see  less  leaf  surface  area  than  will  photons  traveling  at  smaller  angles.  Therefore,  photons 
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propagating  near  the  vertical  will  in  general  experience  a  greater  number  of  leaf  interactions  and  have  a 
smaller  mean  free  path  than  those  traveling  at  greater  angles. 

To  derive  an  expression  for  the  mean  projected  area,  first  suppose  that  the  photon  is  traveling  along  the 
direction  s  which  makes  an  angle  of  6  with  respect  to  the  z-axis,  as  illustrated  in  Figure  27(a).  Because 
we  will  assume  that  the  leaf  angular  distribution  is  uniform  in  the  azimuth  angle  (f) ,  we  can  assume  that 
p  lies  in  the  x-z  plane  and  still  arrive  at  a  general  result.  Therefore  the  photon  propagation  vector  can  be 
expressed  as, 


p  =  xsind  +  ircos#  (6) 

The  leaf  angle  orientation  is  described  by  the  normal  vector  n  ,  which  can  be  written  as, 


n-x sin 0L  cos (f)L+y  sin 0L  sin (f)L+z cos 0, 


(7) 


where  X  ,  y  and  z  are  Cartesian  coordinate  system  unit  vectors.  The  projected  leaf  area  is  then  found  by 
means  of  a  vector  dot  product.  Thus,  the  projected  area  in  the  direction  of  photon  propagation  can  be 
expressed  as, 


AP  =  Aon  *  P  (8) 

where  A0  =  7vd2L  /4  is  the  actual  leaf  area  and  where  we  have  tacitly  assumed  circular  leaves  with  a  mean 

effective  diameter  d L  of  7.5cm.27  Inserting  the  expressions  from  Eq.  (6)  and  Eq.  (7)  into  Eq.  (8),  and 
simplifying,  yields  the  following  result 

A  =  A0  [sin  0,  cos  <f>,  sin  6  +  cos  <f>,  cos  0\  (9) 

The  last  expression  is  the  projected  area  of  a  single  generalized  leaf  described  by  a  specific  set  of  leaf 
orientation  angles.  To  find  the  mean  projected  area  Ap  \0),  we  must  then  average  this  expression  over  all 

possible  orientation  angles.  Note,  though,  that  the  angles  6 L  and  (j),  are  governed  by  their  own 
statistically  independent  probability  density  functions  p0 ,  (0L  )  and  /},„  {(f),  ),  respectively,  and  so  are 
included  in  the  integration  as  follows:28 

co  oo 

Ap  M  =  j  j  Ap  » <t>L  '■  0)PgL  (0L  H,  (10) 

—00—00 

Most  often  the  azimuth  angle  distribution  is  assumed  to  be  uniform  on  the  range  [(),27T ) ,  with  probability 
density  function  (PDF)  then  written  as,  29,30 


PoL  {</>l  ) 


[Ml  )  ~  u{<l>L  -2 7t)\ 

2  n 


(11) 
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where  u(x-x0  MG  if  x  >  x0)or  (Oif  X  <  J0)]  is  the  unit  step  function.31  The  zenith  angle  PDF  p&  is, 
however,  commonly  considered  to  be  cosinusoidal,  being  expressed  as27 


P&L  (oL) 


u(0L)-u 

k-fTI 

V  2JJ 

cos  0L 


(12) 


Notice  that  this  PDF  corresponds  to  leaf  normal  vectors  that  are  primarily  vertical  and  only  rarely 
horizontal.32 

The  expression  for  the  mean  projected  area  can  then  be  found  by  inserting  Eqns.  (9),  (1 1)  and  (12)  into 
Eq.  (10).  Solving  the  integral  and  simplifying  yields  the  following  final  expression  for  the  mean  projected 
area, 


A 


p 


sin#  + 


—  COS0 

4 


(13) 


3.2.I.U.  Leaf  Number  Density 

It  is  commonly  accepted  that  the  volume  distribution  of  leaves  varies  as  a  function  of  vertical  height  z 
within  the  canopy.33'34  Often  the  vertical  profile  of  a  canopy  is  described  in  terms  of  leaf  area  density 

Z(z),  which  is  a  measure  of  the  total  leaf  area  contained  within  a  volume  of  the  canopy.  This  is  a 
convenient  way  to  characterize  the  distribution  as  the  leaf  area  density  is  readily  converted  to  number 
density  by  dividing  by  the  actual  mean  leaf  area. 

The  expression  we  used  to  define  the  leaf  area  density  function  was  empirically  derived  by  other 
researchers  and  is  based  on  several  previously  archived  forest  parameters,  including:  total  canopy  height 
h  ,  maximum  leaf-area  density  Lm ,  and  the  corresponding  canopy  height  zm  at  which  the  leaf  area 
density  takes  on  its  maximum  value.34  In  particular, 


where 


i(d=4 


'h 

V 


z  'j 

m 

h-z  j 


exp 


h  ~  z„ 
h  -  z 


(14) 


j  6  0 <z<zm 

l1/2  Zm^Z^h 


As  an  example,  the  leaf  area  density  for  the  canopy  we  used  in  our  simulation  is  plotted  in  Figure  28, 
where  h  —  18 m ,  Lm  =  0. 6132m-1 ,  and  zm  j h  -  0.85  .  Note  that  the  value  for  Lm  corresponds  to  our 
own  experimental  results,  as  will  be  described  in  Section  4,  while  the  other  values  were  chosen  in 
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accordance  with  maple  tree  forest  parameters  measured  by  others.14  The  leaf  number  density  function 
N(z)  is  then  found  by  dividing  the  leaf  area  density  by  the  actual  leaf  area, 


AJ(Z)  = 


7td~ 


h-  z  ) 


exp 


h-. 


h  —  z 


(15) 


Finally,  the  expression  for  the  mean  free  path  between  scatterers  can  be  determined  by  inserting  Eqns. 
(13)  and  (15)  into  Eq.  (5),  and  simplifying,  to  yield 


D(z,0) 


71 
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sin  0  +  —  cos  0 
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exp 
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3. 2.1. Hi.  Selecting  a  Random  Propagation  Distance  in  an  Inhomogeneous  Medium 

Recall  that  Eq.  (4)  gives  the  probability  density  function  for  the  random  propagation  distance  in  a 
homogeneous  medium.  However,  in  actuality  tree  canopies  are  inhomogeneous,  as  evidenced  by  the  leaf 
number  density  changing  as  a  function  of  canopy  height.  Our  approach  to  account  for  this  fact  was  to 
break  the  canopy  into  fifty  vertically-stacked  horizontal  slices,  so  that  the  homogeneity  of  the  leaf  number 
density  was  approximately  preserved  within  each  slice.  In  order  to  minimize  the  total  number  of  slices 
we  divided  the  canopy  into  non-uniformly  spaced  regions  where  the  absolute  value  of  the  edge-to-edge 
change  in  leaf  area  density  was  common  to  each  region.  As  a  result,  the  canopy  slices  in  regions  where 
the  leaf  area  density  changes  rapidly  are  thinner  than  in  those  regions  where  it  varies  slowly.  The  slices 
were  then  treated  as  separate  canopies  stacked  one  on  top  of  the  other,  and  propagation  through  each  slice 
was  treated  on  a  region  by  region  basis. 

In  any  given  homogeneous  slice  of  the  canopy,  a  random  propagation  distance  d  is  selected  using  the 
mean  free  path  of  that  region.  If  the  propagation  distance  is  less  than  the  distance  necessary  for  the  photon 
to  leave  the  current  region,  then  we  assume  that  there  has  been  a  leaf  interaction.  Otherwise  we  truncate 
the  propagation  at  the  edge  of  the  current  region  and  begin  propagation  again  in  the  next  region.  Note  that 
if  the  random  propagation  distance  would  cause  an  upward  traveling  photon  to  leave  the  top  most  region 
of  the  canopy  we  propagate  the  photon  back  to  the  receiver  plane  and  proceed  as  discussed  earlier. 
Similarly,  if  the  random  propagation  distance  would  cause  a  downward  traveling  photon  to  exit  the 
bottom  most  region,  we  then  terminate  propagation  at  the  ground  plane  and  determine  whether  the  photon 
has  struck  either  the  ground  or  one  of  the  two  targets. 


Figure  9  Representative  Maple  Leaf  Area  Density  as  a  Function  of  Canopy  Height 
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Figure  10  Illustration  of  Region-to-Region  Propagation  in  a  Vertically  Segmented  Tree  Canopy 

For  example,  consider  a  photon  propagating  at  an  angle  9  whose  last  leaf  interaction  occurred 
somewhere  in  the  kth  region,  as  shown  in  Figure  29  (a).  The  vertical  distance  between  the  current  location 
of  the  photon  and  the  boundary  of  the  next  region  is  df  and  the  total  thickness  of  kth  region  is  given  by  tk 

.  Also,  suppose  a  random  value  for  dk  is  chosen  which  places  the  photon  outside  the  kth  region;  that  is, 
dk  •  cos  9  >  df .  We  then  set  the  total  distance  traveled  in  this  jump  to  |  df  /cos  ()\  and  begin  propagating 
again  at  the  edge  of  the  (k+l)th  region.  Next  we  use  the  mean  free  path  for  the  (k+l)th  region  to  select  a 
new  propagation  distance  dk+l ,  which  is  then  compared  to  the  thickness  of  the  (k+l)th  region.  If 

dk+l  •  COS  9  >  tk+1 ,  as  is  illustrated  in  Figure  29  (b),  then  the  cumulative  distance  traveled  is  set  to 
|  df  / cos  9  |  + 1  tk+l  / cos  9  |  and  the  process  is  restarted  at  the  edge  of  the  (k+2)th  region.  However,  if 
dk+l  ■  cos  9  <  tk+ j ,  as  shown  in  Figure  29  (c),  then  the  photon  is  assumed  to  have  experienced  a  leaf 
interaction  in  the  (k+l)th  region,  the  total  distance  traveled  is  set  to  |  df  /cos  9\  +  dk+l  and  the  propagation 
sequence  is  restarted  with  new  scattering  angles. 

3.2.2.  Leaf  Scattering  Angles 

Recall  that  when  a  photon  interacts  with  a  leaf  either  absorption  or  scattering  may  occur.  As  previously 
discussed,  we  determine  whether  or  not  the  photon  is  absorbed  by  selecting  a  value  of  a  random  variable 
uniformly  distributed  between  0  and  1  and  comparing  it  to  the  absorption  coefficient  A  .  In  a  similar 
fashion,  when  scattering  occurs  we  determine  the  form  of  scattering  (i.e.  either  reflection  or  transmission) 
by  also  selecting  a  value  of  a  random  variable  uniformly  distributed  between  0  and  1 .  In  this  case,  though, 
if  the  randomly  selected  value  is  between  0  and  R/{R  +  T^,  where  R  and  T  are  the  reflection  and 
transmission  coefficients  of  the  leaf,  respectively,  then  reflection  is  assumed  to  have  occurred.  Otherwise, 
the  photon  is  assumed  to  be  transmitted.  In  our  previous  work  we  measured  R  ,  T ,  and  A  for  healthy 
sugar  maple  leaves  to  be  0.4861,  0.4841,  and  0.0298,  respectively. 
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Regardless  of  the  form  of  scattering,  the  leaf  scattering  angles  9S  and  </>s  are  generally  assumed  to  be 
statistically  independent  with  separable  PDF’s.35  The  scattering  azimuth  angle  is  most  often  taken  to  be 
uniformly  distributed  on  the  interval  [0,2;r) ,  with  its  probability  density  function  p0 ^  ((/)s  )  written  in  the 

same  form  as  Eq.  (11)  but  with  (j)s  substituted  for  (j)L .  We  have  also  previously  shown  that  in  reflection 

the  zenith  scattering  angle  can  be  modeled  as  the  summation  of  two  probability  distributions:  a  lambertian 
distribution  accurately  models  diffuse  scattering  from  the  leaf,  while  the  specular  reflection  component 
can  be  modeled  by  a  Rayleigh  distribution.  The  normalized  sum  of  the  two  components  is  then  expressed 
as 
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where  Fs  is  the  specular  reflection  fraction  and  0lnc  is  the  leaf  illumination  angle.  In  the  event  that 
transmission  takes  place  the  scattering  distribution  function  is  assumed  to  be  perfectly  Lambertian  and  the 
specular  reflection  fraction  is  set  to  zero.  Thus,  the  expression  for  p&  (9S),  above,  is  valid  for  both 
reflection  and  transmission  and  integrates  to  unit  area  for  each  individually. 

A  sample  bidirectional  scattering  distribution  function  (BSDF)  for  a  sugar  maple  leaf  ( Acer  saccarum), 
fitted  to  data  measured  at  an  incident  zenith  angle  of  1 10°  is  shown  in  Figure  30.  Incident  light,  shown  by 
the  arrow,  hits  the  leaf  at  110°  (or  70°  above  the  horizontal  axis),  while  the  leaf  itself  is  assumed  to  be 
oriented  vertically  so  that  its  surface  normal  vector  n  is  parallel  to  the  0°  axis.  The  BSDF  lobe  on  the  left 
is  then  due  to  reflection  and  the  lobe  on  the  right  is  due  to  transmission.  In  our  previous  work  we 
measured  BSDF’s  for  both  sugar  maple  and  eastern  cottonwood  ( Populus  deltoides)  trees  in  the  local 
Dayton,  Ohio,  area  and  built  extensive  BSDF  models  which  we  have  used  in  our  Monte  Carlo  canopy 
scattering  simulation. 
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Figure  11  Example  BSDF  for  a  Sugar  Maple  Leaf  Illuminated  at  an  Incidence  Zenith  Angle  of  110° 
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3.2.3.  Ground  and  Target  Scattering  Angles 


When  photons  reach  the  ground  plane,  their  location  is  compared  to  the  locations  of  the  targets.  Photons 
which  do  not  hit  either  of  the  targets  are  assumed  to  be  scattered  from  the  ground.  The  azimuth  and  zenith 
ground  scattering  angles  are  again  commonly  assumed  to  be  statistically  independent.  In  particular,  the 

ground  azimuth  scattering  angle  distribution  p0  (<f>c  )  is  assumed  to  be  uniform  over  the  range  [0,2;r) 

and  can  be  expressed  using  Eq.  (11)  after  substituting  (j)G  for  (f)L  .  The  ground  zenith  scattering  angle 

distribution  p@  (0C:  )  is  then  governed  by  a  simple  Lambertian  distribution,  nominally  describing  upward 

scattering  in  the  positive  z-direction,  and  can  be  expressed  using  Eq.  (12).  after  substituting  0G  for  0L  ,36 

Photons  which  do  hit  a  target  are  scattered  according  to  the  probability  density  functions  governing  the 
target  zenith  and  azimuth  scattering  angles.  Once  again,  the  two  random  angles  are  considered  to  be 
statistically  independent.  As  with  ground  scattering,  the  target  azimuth  scattering  angle  distribution 

{(f)T )  is  assumed  to  be  uniformly  distributed  over  the  interval  [0,2;r )  and  the  target  zenith  scattering 
angle  distribution  p&  {0T  )  is  taken  to  be  Lambertian.  The  PDFs  p,b  and  pG)/  can  then  be  found  by 
substituting  (j)T  for  (f)L  ,  and  0T  for  0L  in  Eqns.  (11)  and  (12),  respectively. 

Recall,  however,  that  both  targets  are  assumed  to  be  tilted  so  that  their  surface  normals  are  pointed  back 
in  the  direction  of  the  receiver.  As  a  result,  the  target  scattering  angles  must  be  transformed  into  the 
canopy  coordinate  system  as  previously  discussed.  We  have  done  this  for  two  reasons,  the  first  being  that 
in  our  simulation  we  increase  the  probability  of  receiving  photons  back  from  the  targets  if  we  nominally 
tilt  them  back  towards  the  receiver.  And  secondly,  in  real  world  situations  targets  will  typically  have  one 
or  more  facets  that  do  in  fact  face  back  toward  the  receiver. 

3.2.4.  Test  Canopy  Parameters 

We  wished  to  verify  the  results  of  our  canopy  propagation  model  with  real  data  collected  from  within  a 
nearby  grove  of  trees.  Therefore,  in  order  to  achieve  the  most  reproducible  results,  we  used  real 
parameters  measured  for  this  tree  grove  in  our  simulation.  These  include  the  canopy  dimensions,  the 
illumination  angle  geometry,  and  the  leaf  area  density. 

We  illuminated  our  tree  grove  using  a  rudimentary  direct  detection  ladar  system  located  on  the  1 1th  floor 
of  a  tower  building  at  Wright-Patterson  Air  Force  Base,  OH.  Our  ladar  system  employs  a  Raman  shifted 
Nd:YAG  Coherent  Infinity  pulsed  laser  that  has  been  Raman  shifted  to  operate  at  a  wavelength  of  1560 
nm  and  produces  pulses  of  ~1  ns  duration  at  a  rate  of  10  Hz.  We  frequency  doubled  the  operating 
wavelength  with  a  LiNbCE  crystal  to  780  nm  and  added  telescoping  optics  to  control  the  beam  diameter 
down  range.  The  detector  is  a  Princeton  Instruments  Pi-Max2  intensified  CCD  camera  which  is 
synchronized  with  the  laser  for  range  gating  applications.  Using  a  differential  GPS  unit,  we  logged  the 
locations  of  the  primary  target  site  and  the  ladar  system  located  in  the  tower.  We  then  determined  the 
angle  made  between  the  two  locations,  as  well  as  the  distance  separating  them.  We  found  that  the 
illumination  angle  was  approximately  80°  and  the  horizontal  distance  between  the  tower  and  target  was 
255  meters. 

We  then  measured  the  dimensions  of  the  tree  grove  by  walking  around  the  perimeter  with  a  handheld 
GPS  device  and  periodically  logging  waypoints  around  the  edge  of  the  tree  line.  A  best  fit  ellipse  was 
then  drawn  around  the  perimeter,  with  major  and  minor  axis  dimensions  found  to  be  358m  and  228m, 
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respectively.  Conveniently  for  us,  both  the  target  site  and  the  tower  were  found  to  lie  along  the  major  axis 
of  this  ellipse,  allowing  us  to  easily  incorporate  the  actual  canopy  dimensions  into  our  simulation. 

Next  we  indirectly  computed  the  maximum  leaf  area  density  Lm  by  first  measuring  the  gap  fraction  of 

our  test  forest  at  80°,  and  by  then  finding  the  maximum  leaf  area  density  value  which  would  yield  the 
same  statistical  probability  of  a  photon  passing  through  the  canopy  without  interaction.  To  measure  the 
gap  fraction,  we  began  by  capturing  sixty  hemispheric  images,  similar  to  the  one  shown  in  Figure  1,  from 
a  5x12  point  grid  of  locations  under  the  canopy.  We  then  applied  a  threshold  to  the  images  so  that  pixels 
containing  leaves  and  branches  were  given  a  value  of  unity  and  those  without  were  given  a  value  of  zero. 
Each  image  was  then  broken  into  nine  equal  area  annuli  whose  bin  centers  have  been  mapped  to  non- 
uniformly  spaced  zenith  angles.37  The  number  of  zero  valued  pixels  was  counted  and  then  divided  by  the 
total  number  of  pixels  for  each  annulus  in  order  to  determine  the  gap  fraction  as  a  function  of  zenith 
angle.  From  the  60  individual  images,  mean  gap  fraction  values  were  then  calculated  and  plotted  at  the 
center  of  each  angular  bin,  as  shown  in  Figure  31.  The  interpolated  gap  fraction  value  corresponding  to  80 
degrees  was  found  to  be  1 1.7e-3. 

Gap  fraction,  typically  defined  as  the  percentage  of  canopy  area  not  covered  by  leaves,  also  describes  the 
probability  of  encountering  a  gap  at  a  specific  zenith  angle.  That  is,  if  a  photon  is  launched  into  the 
simulation  at  a  certain  zenith  angle,  the  probability  that  it  will  reach  the  ground  without  interaction  is 
equal  to  the  gap  fraction  of  the  canopy  at  that  angle.  Because  we  modeled  propagation  through  the  canopy 
as  a  statistical  process,  we  can  therefore  also  calculate  the  probability  of  a  photon  passing  unscattered 
through  the  simulated  canopy  as  a  function  of  leaf  area  density. 

Recall  that  the  leaf  number  density  varies  as  a  function  of  canopy  height  and  that  for  our  model  we  have 
broken  the  canopy  into  50  regions  of  constant  density.  Also  recall  from  Eq.  (4)  that  we  model  the  photon 
propagation  distance  between  scattering  events  using  a  negative  exponential  probability  density  function. 
Therefore,  the  probability  that  a  photon  traveling  at  zenith  angle  9  experiences  a  leaf  interaction  within 
the  kth  region,  (which  has  mean  free  path  Dk  and  upper  and  lower  bounds  zkl  and  zk  ,  respectively,  as 

shown  in  Figure  29  (c)  can  be  found  by  first  integrating  this  PDF  over  the  region  boundaries  and  then 
dividing  by  a  normalization  factor  according  to 


where  ak  =  zk  /cos  0  .  We  normalize  the  above  slant  path  integral  in  order  to  account  for  the  fact  that  if 

the  photon  experiences  its  first  leaf  interaction  within  the  kth  region,  then  it  must  have  already  passed 
through  the  previous  k- 1  regions  without  scattering. 

The  probability  that  there  will  be  no  leaf  interaction  within  the  kth  region  is  then  found  by  subtracting  Eq. 
(18)  from  unity.  Correspondingly,  the  probability  of  not  having  an  interaction  within  the  entire  canopy 
can  be  found  by  evaluating  the  following  product  relationship 
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Figure  12  Mean  Gap  Fraction  as  a  Function  of  Zenith  Angle  for  our  Nearby  Gove  of  Sugar  Maple  Trees 


PjO  <z  <h\=Y | 


k=\ 


(19) 


We  evaluated  Eq.  (19)  for  an  illumination  angle  of  80°  as  we  varied  the  maximum  leaf  area  density. 
(Recall  that  the  mean  free  path  is  a  function  of  Lm  ,  as  given  by  Eq.  (16).  The  result  is  shown  in  Figure 
32,  and  the  equation  of  the  best  fit  curve  was  found  to  be 


i5{o<x</i}=exp(-  9.0490  ■  Lm)  (20) 

From  this  result,  we  found  that  a  maximum  leaf  area  density  of  Lm  =  0.6 1 32  m2/  m  would  yield  the 

same  probability  of  an  unscattered  photon  reaching  the  ground  as  was  determined  from  the  gap  fraction 
measurements.  This  value  is  shown  by  the  arrow  in  Figure  32. 


3.3  Experimental  Verification 

In  order  to  validate  our  model,  we  collected  waveform  data  from  a  nearby  stand  of  sugar  maple  trees 
during  the  Summer  of  2008.  We  performed  waveform  measurements  by  illuminating  the  canopy  at 
approximately  a  10°  angle  above  the  horizon  (i.e.,  80°  zenith  angle)  and  then  sampling  the  beam  footprint 
on  the  ground  with  a  grid  of  avalanche  photodiodes  (APDs).  We  sent  two  beams  from  the  tower,  one  at 
780nm  and  the  other  at  1560nm.  An  image  of  the  test  canopy  as  viewed  from  the  tower  is  shown  in 
Figure  38,  where  the  beam  was  aimed  at  a  target  placed  near  the  1C  location. 
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Figure  13  Image  of  the  Test  Canopy  as  Viewed  from  the  Tower 


We  established  a  reliable  trigger  signal  by  placing  aim  diameter,  high  reflectance  spectralon  target  in  a 
clearing  so  that  it  was  clearly  visible  from  the  tower.  Thus  we  illuminated  with  the  1560  nm  beam.  An 
amplified  PIN  diode  was  then  stationed  0.5  m  in  front  of  the  spectralon  target  and  aimed  in  such  a  way 
that  the  entire  target  fell  within  the  field  of  view  of  the  detector.  We  then  relayed  the  trigger  signal 
through  a  500  m  co-axial  cable  to  an  oscilloscope  stationed  under  the  canopy.  Keeping  the  trigger  signal 
at  a  constant  location  allowed  synchronization  of  the  temporal  delays  between  waveforms  measured  at 
different  grid  point  locations. 

We  then  illuminated  the  tree  grove  with  the  780nm  beam  which,  due  to  diffraction,  had  a  width  of  about 
10m  at  the  top  of  the  canopy.  This  beam  was  directed  towards  a  reflexite  (primary)  target  located  near  the 
center  of  the  grove.  Note  that  although  our  model  uses  maple  leaf  BSDF’s  measured  at  1064nm,  the 
scattering  functions  of  these  leaves  are  expected  to  be  very  similar  under  780nm  illumination.  We  then 
cleared  an  88m  path  from  the  primary  target  toward  the  tower,  with  twelve  additional  16m  paths  cleared 
orthogonal  to  the  beam  propagation  axis  and  spaced  in  uniform  8m  intervals.  Along  these  paths  we  then 
established  five  uniformly  spaced  measurement  locations  to  create  a  12x5  measurement  grid,  as  seen  in 
Figure  38.  An  image  of  the  cleared  path  beneath  the  canopy  can  be  seen  in  Figure  39.  We  observed  that 
the  beam  appeared  to  enter  the  canopy  roughly  above  the  eleventh  transverse  path,  approximately  80m 
from  the  primary  target. 


Figure  14  12x5  Measurement  Grid  Beneath  the  Canopy 
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Figure  15  Image  of  the  Path  Connecting  the  Entrance  Location  to  Target  as  Seen  from  the  Entrance  Location 


We  then  measured  the  waveforms  at  each  grid  location  with  a  Hamamatsu  silicon  avalanche  photodiode, 
seen  in  Figure  41,  which  was  mounted  on  a  tripod,  aimed  directly  upwards,  and  moved  from  point  to 
point  throughout  the  test.  This  detector  had  a  field  of  view  of  approximately  45°.  We  then  relayed  the 
detected  signals  to  the  oscilloscope  along  with  the  trigger  signal.  At  each  grid  point  we  collected  256 
pulses,  which  we  averaged  into  a  single  mean  wavefonn  and  stored  for  later  processing. 
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4.  Results  and  Discussions 


4.1  Leaf  BSDF  Measurement  Results 

Polar  plots  of  the  BSDF  data  collected  from  maple  leaves  in  the  manner  described  above  are  provided  in 
Figure  5.  The  laser  beam  (depicted  by  the  arrows)  was  incident  upon  the  leaves  at  0,  angles  of  0°,  10°, 
20°,  30°,  40°,  50°,  60°,  70°,  and  78°.  A  separate  sub-figure  was  made  for  each  illumination  angle.  Each 
sub-plot  has  been  normalized  such  that  the  sum  of  the  areas  of  the  transmission  and  reflection  lobes  is 
unity.  Note  that  the  scaling  of  the  radial  axis  spans  a  range  from  0  to  0.01  for  each  subplot,  with  exception 
of  78°  illumination  which  goes  from  0  to  0.025. 

The  shape  of  the  BRDF  and  the  BTDF  at  normal  incidence  and  low  illumination  angles  appear  to  be 
largely  Lambertian  in  nature.  The  radiated  energy  is  located  mostly  along  the  leaf  surface  normal, 
decreasing  cosinusoidally  with  the  rotation  of  the  detector  angle.  In  addition,  the  shape  of  the  BTDF 
remains  relatively  constant  regardless  of  incident  angle.  However,  as  the  incident  angle  increases  past  50°, 
the  specular  component  of  the  reflection  increases  and  the  transmission  decreases.  The  glint  protrudes 
only  slightly  from  the  diffuse  component  at  50°  and  becomes  more  pronounced  as  the  incidence  angle 
increases.  At  78°,  the  specular  component  dominates  the  diffuse  reflection.  Similar  trends  are  seen  in  the 
cottonwood  BSDF  data  shown  in  Figure  6. 

4.1.1.  Calculation  of  Absorption  Coefficient 

Before  investigating  BRDF  and  BTDF  features  in  more  detail,  it  is  important  to  determine  the  absoiption 
coefficient  Al(Oi)  as  a  function  of  illumination  angle.  The  method  involves  measuring  the  BRDF  of  a 
target  with  a  known  reflection  coefficient  Rs  and  comparing  its  area  to  that  of  the  BRDF  and  BTDF  of  a 
leaf.  A  60%  reflective  Spectralon®  disc  (i.e.,  Rs=0.60  for  all  0j)  was  used  as  the  standard  of  measure  to 
which  the  scattering  distributions  of  the  leaves  were  compared.  The  normal  illumination  BRDF  of  the 
Spectralon®  disc  is  plotted,  for  example,  along  with  the  maple  leaf  BSDF  for  normal  illumination  in 
Figure  7.  The  area  of  each  function  is  then  calculated  by  integrating  over  the  span  of  detector  angles.  That 
is,  the  area  of  the  maple  leaf  BRDF  is  given  by 


+270 

ArM  =  {  rL  (&,  ,  dD  )ddD 
+90 


(21) 


where  rL(0i ,  0D )  is  the  BRDF  data  measured  at  illumination  angle  0/  and  detector  angle  0D.  Similarly,  the 
area  of  the  BTDF  is  given  by 

+90 

Ar(6>i)=  \t/(0n0D)d0D  (22) 

-90 


where  ti/Oj ,  0D)  is  the  BTDF. 
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Figure  16  Measured  BSDF  Data  of  Common  Maple  Leaves  for  Illumination  Angles  6j  of  (a)  0°  (b)  10°  (c)  20° 

(d)  30°  (e)  40°  (f)  50°  (g)  60°  (h)  70°  and  (i)  78° 
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Figure  17  Measured  BSDF  Data  of  Cottonwood  Leaves  for  Illumination  Angles  6j  of  (a)  0°  (b)  10°  (c)  20°  (d) 

30°  (e)  40°  (f)  50°  (g)  60°  (h)  70°  and  (i)  78° 
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The  values  of  the  reflection  RL(0i )  and  transmission  Tr  ( 0/,  coefficients  are  directly  related  to  the  areas 
under  the  BRDF  and  BTDF  curves.  In  particular,  the  ratio  of  the  reflection  coefficients  for  maple  leaves 
and  the  Spectralon®  disc  is  equal  to  the  ratio  of  these  areas  according  to  the  relationship: 


Rl(°i)  _  Ar(°i)  (23) 

0.60  As 

where  As  is  the  area  of  the  spectralon  BRDF  at  normal  illumination.  A  similar  expression  can  be  written 
for  the  transmission  coefficient: 


tM  aM 

0.60  As 

The  absorption  coefficient  can  then  be  found  through  the  expression 


(24) 


(25) 


Table  1  and  Table  2  present  values  for  the  reflection,  transmission,  and  absorption  coefficients  for  maple 
and  cottonwood  leaves  as  a  function  of  incident  angle.  Notice  that  the  transmission  coefficients  decrease 
and  the  reflection  coefficients  increase  with  increasing  illumination  angle.  This  trend  is  also  seen  in  the 
polar  plots  of  the  leaf  BSDFs  in  Figures  5  and  6  where  the  size  of  the  transmission  lobe  is  seen  to  grow 
smaller  with  increasing  illumination  angle.  Another  trend  evident  in  Tables  1  and  2  is  the  growth  of 
absorption  with  illumination  angle. 


Figure  18  Comparison  of  the  Sugar  Maple  Leaf  BSDF  and  Spectralon®  BRDF  for  Normal  Illumination 
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Figure  19  Absorption  Coefficient  for  Maple  Leaves  as  a  Function  of  Illumination  Angle 


The  relationship  between  illumination  angle  and  the  absoiption  coefficient  for  maple  leaves  is  illustrated 
in  Figure  8.  A  second  order  polynomial  was  fit  to  the  data  in  order  to  allow  the  absoiption  coefficient  for 
any  illumination  angle  to  be  calculated.  The  following  two  equations,  then,  describe  least  square  fit 
regression  equations  for  the  absoiption  coefficients  of  maple  and  cottonwood,  respectively, 


A(A)  =  -4.5156x10"6  +  1.0373xl0"3  x#7  +25.415xl0"3  (Zb) 

A (^/ ) ~  1  •  1885 x  1  O'6  x 6^  - 67.238 x  10'6  x0;  +32.7xl0  3  (27) 


If  desired,  the  reflection  and  transmission  coefficients  can  be  determined  for  any  illumination  angle  by 
integrating  the  BRDF  and  BTDF  models  we  will  develop  in  the  following  sections. 


Table  1  Sugar  Maple  Leaf  Reflection,  Transmission  and  Absorption  Coefficients  as  a  Function  of 

Illumination  Angle 


Rl 

tl 

al 

0 

0.4861 

0.4841 

0.0298 

10 

0.4879 

0.4784 

0.0337 

20 

0.4890 

0.4727 

0.0383 

30 

0.4884 

0.4661 

0.0456 

40 

0.4900 

0.4406 

0.0694 

50 

0.4984 

0.4300 

0.0716 

60 

0.5000 

0.4327 

0.0673 

70 

0.5009 

0.4235 

0.0756 

78 

0.5031 

0.4188 

0.0781 

Table  2  Eastern  Cottonwood  Leaf  Reflection,  Transmission  and  Absorption  Coefficients  as  a  Function  of 

Illumination  Angle 


Rl 

TL 

al 

0 

0.5414 

0.4263 

0.0323 

10 

0.5451 

0.4227 

0.0322 

20 

0.5500 

0.4178 

0.0322 

30 

0.5519 

0.4160 

0.0321 

40 

0.5548 

0.4131 

0.0321 

50 

0.5596 

0.4084 

0.0320 

60 

0.5457 

0.4221 

0.0322 

70 

0.5836 

0.3847 

0.0317 

78 

0.5801 

0.3834 

0.0365 
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4.1.2.  Modeling  Specular  Reflection 

In  Figure  20  (a),  the  maple  leaf  BRDF  and  BTDF  for  normally  incident  light  are  fitted  with  a  cosine  curve 
normalized  to  the  integrated  area  of  the  measured  data.  The  measured  data  is  depicted  by  the  circles, 
while  the  fitted  curve  is  given  by  the  solid  line.  The  quality  of  the  fit  suggests  that  both  the  BRDF  and 
BTDF  are  accurately  modeled  by  a  Lambertian  distribution  function  for  normal  incidence.  Similar  results 
were  seen  for  the  other  illumination  angles  less  than  50°.  However,  as  the  illumination  angle  increases 
beyond  50°,  other  features  become  apparent.  For  example,  plotted  in  Figure  20  (b)  are  the  70°  illumination 
angle  BRDF  and  BTDF  data  fitted  with  normalized  area  cosine  curves.  While  the  transmission  data 
remains  nearly  Lambertian,  the  reflection  data  exhibits  other  features  that  will  be  addressed  in  this 
section. 

Because  of  the  irregularities  in  the  shape  of  the  reflection  curves,  it  is  difficult  to  model  the  BRDF  data 
with  a  simple  distribution.  Also,  the  specular  reflection  peak  at  high  incident  angles  makes  it  difficult  to 
accurately  replicate  the  curves  by  simply  using  high  order  polynomials.  One  approach  to  modeling  the 
BRDF  is  to  separate  the  specular  and  diffuse  components.  Accurate  models  individually  representing  the 
diffuse  and  specular  profiles  can  then  be  found  separately  and  later  added  together  to  produce  the  full 
BRDF  model. 

Under  the  assumption  that  the  specular  component  is  negligible  for  detector  angles  less  than  40°,  the 
BRDF  can  be  broken  into  two  regions:  one  corresponding  only  to  diffuse  reflection  ( 0D  <40 ),  and  the 
other  containing  both  specular  and  diffuse  elements  {0D  >  40).  Separation  of  the  two  components  in  the 
latter  region  is  then  performed  by  fitting  a  polynomial  curve  to  the  known  diffuse  data  and  interpolating 
values  within  this  region.  Subtracting  the  interpolated  diffuse  data  from  the  overlapping  region  leaves 
only  the  specular  reflection  component. 


(a) 


(b) 


Figure  20  Lambertian  Fit  to  the  BRDF  (left)  and  BTDF  (right)  for  (a)  Normal  Illumination  and  (b) 

Illumination  at  70° 
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Figure  21  Separation  of  the  Sugar  Maple  Leaf  Diffuse  and  Specular 
Reflection  Components  for  70°  Illumination 

This  process  is  illustrated  in  Figure  21  for  maple  leaf  BRDF  data  taken  at  70°  illumination  angle.  The 
non-separated  BRDF  data  is  depicted  by  the  stars  and  circles,  which  in  turn  represent  the  diffuse  and 
overlapping  regions,  respectively.  A  fourth  order  polynomial  is  fit  to  the  diffuse  data  and  represented  by 
the  solid  line.  The  specular  reflection  component  is  then  found  by  subtracting  the  diffuse  reflection  fit 
from  the  measured  data.  This  is  shown  for  the  50,  60,  70,  and  78  degree  illumination  angles  in  Figure  22. 

As  also  shown  in  Figure  22,  we  found  that  the  specular  data  is  best  fit  with  a  normalized,  reversed 
Rayleigh  distribution  of  the  form, 


rM=FM)S^fexp 


(90  ~  0D  )2 
2(90  -  dP  )2 


(28) 


where  Fs  {(), )  is  the  fractional  specular  reflection,  0D  is  the  detector  angle,  and  0P  is  the  angle  at  which 

the  function  is  a  maximum.  Interestingly,  we  found  that  a  6P  value  of  75°  produced  the  best  fit  for  all 
incidence  angles. 

The  area  of  the  reversed  Rayleigh  distribution  must  be  normalized  to  the  fractional  area  of  the  specular 
reflection  peak  before  reconstructing  the  data  with  this  model.  Dividing  the  area  of  the  specular  reflection 
peak,  Ars  (0, ) ,  by  the  total  area  of  the  reflection  lobe,  A  K  (0I ) ,  for  a  given  incident  angle  0,  gives  the 
fractional  specular  reflection  according  to 


=  (29) 

As  is  seen  in  Figure  22,  the  fraction  of  specular  reflection  is  dependent  on  the  illumination  angle.  For 
illumination  angles  less  than  50°,  specular  reflection  is  negligible  and  Fs{dI  <  50" j  =  0 .  Moreover,  at 

90°  incidence  it  is  assumed  that  all  reflection  is  specular,  yielding  Fs  (90°)  =  1 .  Using  these  bounds,  the 

fractional  specular  reflection  was  fit  with  a  third  order  polynomial  as  shown  in  Figure  23  for  maple 
leaves.  The  ratio  of  specular  to  total  reflection  for  any  illumination  angle  0f  (in  degrees)  can  then  be 
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written  empirically  as  a  piecewise  continuous  function  according  to  the  following  equations  provided  for 
Sugar  Maple  and  Eastern  Cottonwood  leaves,  respectively, 


f  0  ,6,  <50  (30) 

[12.7  xlO'6  x  Q]  - 1.7  x  10'3  x  6]  +  77.2  x  10'3  x  6l  -1.19  ,01  >50 

J  0  ,0,  <50 

[- 1.65  x  10'6  x0f  +922  xlO'6  *0]  - 192  xlO'3  x^  +1.89  ,0,  >50 
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Figure  22  Normalized,  Reverse  Rayleigh  Fits  of  the  Specular  Reflection  Data  for  Illumination  Angles  of  (a) 

50°,  (b)  60°,  (c)  70°,  and  (d)  78°. 
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Figure  23  Fractional  Specular  Reflection  for  Maple  Leaves  as  a  Function  of  Incident  Angle 


4.1.3.  Modeling  Transmission  and  Diffuse  Reflection 

Both  the  diffuse  reflection  and  the  transmission  distribution  functions  are  simply  fit  with  polynomials 
such  that  when  added  to  the  Rayleigh  model  for  specular  reflection,  the  complete  BSDF  is  reconstructed. 
Figure  24  shows  the  fitted  data  at  an  illumination  angle  of  70°  for  maple  leaves.  The  transmission  data  is 
fit  with  a  second  order  polynomial,  while  a  fourth  order  polynomial  is  used  for  the  diffuse  reflection.  A 
higher  order  polynomial  is  needed  for  diffuse  reflection  because  the  structure  of  the  distribution  becomes 
somewhat  more  complex  at  higher  incident  angles.  When  the  specular  reflection  component  is  then  added 
to  the  corresponding  polynomials  modeling  diffuse  reflection,  the  complete  BRDF  is  reconstructed,  as 
shown  for  maple  leaves  in  Figure  25. 


Figure  24  Fitted  Sugar  Maple  Leaf  BSDF  Curves  for  Illumination  at  70° 
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Figure  25  Sugar  Maple  Leaf  BRDF  Models  Created  by  Adding  the  Diffuse  and  Specular  Components  at  (a) 

50°,  (b)  60°,  (c)  70°,  and  (d)  78°  Illumination. 


4.1.4.  Leaf  Data  Interpolation 

The  ability  to  accurately  estimate  the  BSDF  for  any  illumination  angle  will  be  a  valuable  resource.  Thus 
far,  it  has  been  shown  that  scattering  data  at  the  measured  illumination  angles  can  be  reconstructed  using 
simple  polynomial  fits  and,  for  high  illumination  angle  BRDF’s,  a  reversed  Rayleigh  distribution. 
However,  filling  in  the  gaps  for  intermediate  illumination  angles  requires  creating  an  additional  fit  to  the 
modeled  data.  As  discussed  in  the  previous  two  sections,  for  a  given  incidence  angle,  the  polynomial 
equations  used  to  describe  the  BTDF  and  diffuse  BRDF  as  a  function  of  both  detector  0D  and 
illumination  Ol  angles  are  given,  respectively,  by: 


t{0D  ,0,)  =  pu0l  +  p2tdD  +  p3t  (  } 

rd(0D>0l)=  Pl,K  +  P2r°l  +  Pir°l  +  Pa^D  +P5r  (33) 

where  the  polynomial  coefficients,  p  ,  are  dependent  on  the  illumination  angle.  With  knowledge  of  these 
coefficients,  the  BTDF  and  diffuse  BRDF  data  can  be  interpolated  at  intermediate  detector  angles  for  a 
previously  examined  illumination  angle.  However,  in  order  to  interpolate  BTDF  and  diffuse  BRDF  data 
at  intermediate  illumination  angles  it  is  necessary  to  examine  the  relationship  between  9j  and  the  value  of 

each  of  the  p-coefficients.  Fitting  p  as  a  function  of  0,  will  allow  the  interpolation  of  p  at  intermediate 
illumination  angles. 

In  Figure  26,  the  p-coefficients  of  the  maple  leaf  BTDF  data  are  shown  fit  with  a  fourth  order  polynomial. 
In  general,  the  fourth  order  polynomial  fits  to  the  p-coefficients  for  both  BTDF  and  diffuse  BRDF  data  is 
defined  by  a  set  of  five  q-coefficients  according  to 
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Pix^Pl)  Ch@ l  +  4 @1  “*”^5 


(34) 


where,  for  transmission,  i  =  1,2,3  and  x  =  t ,  and  for  reflection,  i  =  1,2, 3, 4, 5  and  x  =  r .  Knowledge  of 
the  q-coefficient  values  then  allows  one  to  calculate  the  p-coefficients  for  any  illumination  angle  through 
Eqs.  (2.12)  and  (2.13).  The  q-coefficients  calculated  for  the  transmission  and  diffuse  reflection  Sugar 
Maple  leaf  data  are  shown  in  Table  3  and  Table  4,  respectively.  Similarly,  values  for  the  Eastern 
Cottonwood  q-coefficients  are  provided  in  Table  5  and  Table  6. 


Using  this  method  of  data  interpolation  we  have  examined  the  BSDF  estimates  at  each  of  the  angle 
permutations  used  during  the  data  acquisition  procedure.  The  RMS  errors  found  between  the  originally 
measured  data  and  the  estimates  generated  through  this  interpolation  method  were  found  to  be 
approximately  2.5%  and  1.0%  for  the  BRDF  and  BTDF  of  Sugar  Maple  leaves,  respectively.  The 
corresponding  RMS  errors  calculated  for  Eastern  Cottonwood  leaves  were  approximately  3.7%  and  1.2%. 


(a)  (b)  (c) 

Figure  26  Polynomial  Fits  of  the  (a)  0th,  (b)  1st,  and  (c)  2nd  Order  p-Coeflicients  Describing  the  BTDF  Data  for 

Sugar  Maple  Leaves. 


Table  3  Sugar  Maple  leaf  BTDF  q-Coeflicients 


Qlt 

q2t 

q3, 

q« 

q5t 

Pit 

3.0388e-12 

-463.75e-12 

23.428e-09 

-300.83e-09 

-3.4701e-06 

P2t 

34.837e-12 

-7. 1221e-09 

517.86e-09 

-17.93e-06 

135.77e-06 

P3t 

-10.371e-09 

1.6844e-06 

-89.529e-06 

1.0167e-03 

246.48e-03 

Table  4  Sugar  Maple  leaf  BRDF  q-Coeflicients 


qn 

q  2r 

q3r 

q4r 

q5r 

Plr 

1.2954e-15 

-180.66e-15 

6.8481e-12 

-132.35e-12 

120. 1 6e- 12 

P2r 

25.825e-15 

-5.0125e-12 

285.57e-12 

-3.643  le-09 

32.217e-09 

P3r 

-10.817e-12 

1.5674e-09 

-58.687e-09 

1.0776e-06 

-34.53  le-06 

P4r 

-123.8e-12 

23.992e-09 

-1.2636e-06 

8.0717e-06 

-164.4e-06 

P5r 

-6.2703e-09 

369.59e-09 

-12.501e-06 

-257.65e-06 

245.19e-03 
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Table  5  Eastern  Cottonwood  leaf  BTDF  q-Coefllcients 


Qn 

421 

q3t 

44t 

q5t 

Pit 

2.5913e-12 

-391.07e-12 

17.815e-09 

-191.88e-09 

-34.326e-06 

P2t 

4.4203e-12 

-4. 1417e-09 

325.17e-09 

-4.7745e-06 

-168.86e-06 

P3t 

-9.2296e-09 

1.499e-06 

-75.217e-06 

738.91e-06 

245.53e-03 

Table  6  Eastern  Cottonwood  leaf  BRDF  q-Coefficients 


4lr 

42r 

q3r 

q4r 

q5r 

Pn 

2.1389e-15 

-272.67e-15 

8.6857e-12 

-11 1.72e-12 

1.6618e-09 

P2r 

-59.223e-15 

7.1463e-12 

-208.54e-12 

1.3098e-09 

37.969e-09 

P3r 

-18.62e-12 

2.5167e-09 

-85.345e-09 

1.098e-06 

-47.705e-06 

P4r 

479.27e-12 

-62.848e-09 

2.2671e-06 

-23.813e-06 

-355.59e-06 

P5r 

5.4467e-09 

-1.4462e-06 

62.194e-06 

-993.25e-06 

261.58e-03 

4.1.5.  Procedure  for  Constructing  the  BSDF 

BSDF  values  at  any  set  of  angles  (  0, ,  6 n  1  can  be  accurately  estimated  for  both  maple  and  cottonwood 

leaves  from  the  information  provided  herein.  Using  the  following  procedure,  surface  fits  for  both  the 
BRDF  and  BTDF  can  be  interpolated  for  illumination  angles  spanning  a  range  from  0°  to  78°.  These  are 
plotted  for  maple  leaves,  for  example,  in  Figure  27  and  Figure  28,  respectively. 

1)  Using  Eqs.  (26)  or  (27),  calculate  the  absorption  coefficient,  AL  (O, ) 

2)  Using  Eqs.  (32)  and  (34),  calculate  the  BTDF,  t{dD,0I).  The  coefficients  required  for  these 
equations  are  found  in  Table  3  and  Table  5  for  maple  and  cottonwood  leaves,  respectively. 

3)  Using  Eqs.  (33)  and  (34),  calculate  the  diffuse  BRDF,  rd  (0D  ,0j).  The  coefficients  required  for 
these  equations  are  found  in  Table  4  and  Table  6  for  maple  and  cottonwood  leaves,  respectively. 

4)  Generate  the  specular  BRDF,  r  (0D  ,0,),  using  Eq.  (28).  Normalize  the  function  to  the  fractional 

specular  reflection  calculated  from  either  Eq.  (30)  or  (3 1)  for  maple  and  cottonwood  leaves, 
respectively. 

5)  Construct  the  complete  BRDF  by  adding  together  the  specular  and  diffuse  components.  That  is, 

=  rd(0D>0i)  +  r,{0D>0i)  (35) 
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Scattering  Angle  -100  0  Incident  Angle 
Figure  27  Two  Dimensional  BRDF  Surface  Fit  for  Sugar  Maple  Leaves. 


Scattering  Angle  -100  0  Incident  Angle 
Figure  28  Two  Dimensional  BTDF  Surface  Fit  for  Sugar  Maple  Leaves 


Note  that  the  BRDF  and  BTDF  have  been  normalized  in  such  a  way  that  the  sum  of  their  areas  plus  the 
absorption  coefficient  calculated  in  step  (1)  is  equal  to  unity.  That  is, 

90  270 

1  =  j  t{0D ,  6,  )d 6 n  +  J  rL  {0D ,  6I  )d 0D  +  A,  (0, )  (36) 

-90  90 

The  above  steps  present  a  stand-alone  method  for  generating  the  BSDF  of  maple  and  cottonwood  leaves 
for  any  illumination  angle.  This  method  is  appropriate  for  use  in  remote  sensing  models  where  probability 
density  functions  describing  the  scattering  by  individual  leaves  must  be  considered  (e.g.,  in  Monte  Carlo 
simulations).  Such  models  typically  require  examination  of  the  scattering  PDFs  in  both  the  zenith  and 
azimuth  angles.  While  our  models  have  been  developed  based  only  upon  in-plane  measurements,  they 
may  be  extended  to  allow  out  of  plane  estimates  by  assuming  azimuthal  symmetry  in  both  the  diffuse 
BRDF  and  the  BTDF.  We  address  the  asymmetry  of  the  specular  component  by  forcing  the  azimuth  angle 
to  be  equal  to  the  angle  of  incidence  in  the  case  of  specular  reflection.  We  furthermore  believe  our  models 
to  be  largely  independent  of  polarization,  though  this  will  require  further  investigation  to  verify. 

4.2  Monte  Carlo  Simulation  Results 

Once  the  canopy  parameters  from  the  test  forest  were  incorporated  into  our  model  we  ran  the  Monte 
Carlo  simulation  and  collected  spatial,  temporal,  and  angular  data  in  both  the  ground  and  receiver  planes. 
Over  the  period  of  several  weeks,  we  simulated  the  launching  of  ten  billion  photons  into  the  canopy  at  an 
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Arrival  Time  [ns] 


Figure  29  Characteristic  Temporal  Waveform  from  those  Simulated  Photons  which  Strike  the  Canopy  Floor 


illumination  angle  of  80°.  This  ensured  that  a  sufficient  number  of  photons  returned  back  to  the  detector 
plane,  allowing  us  to  create  smooth  probability  density  functions.  Data  for  all  photons  that  struck  the 
ground  place  were  considered  for  the  ground  plane  analysis,  while  only  those  photons  which  returned  to 
the  detector,  whether  they  hit  the  ground  or  not,  were  considered  in  the  receiver  plane  data  analysis. 

We  first  divided  the  virtual  ground  plane  into  a  grid  of  100x100  rectangular  bins  (i.e.,  the  range  and 
cross-range  bin  dimensions  were  3.58m  and  2.28m,  respectively)  and  created  temporal  waveforms,  such 
as  the  one  shown  in  Figure  30,  from  the  simulated  photons  arriving  within  a  45°  field  of  view.  This  field 
of  view  was  selected  to  match  that  of  the  detector  used  in  our  experiment,  as  we  will  discuss  later  in  this 
section.  These  waveforms  are  generally  characterized  by  a  shaip  peak  due  to  first  surface  scattered 
photons,  followed  by  a  slowly  decaying  tail  due  to  the  arrival  of  multiply  scattered  photons.  We  then 
calculated  the  root  mean  square  (RMS)  pulse  width  of  each  waveform  and  also  calculated  the  integrated 
number  of  photons  falling  within  each  bin.  The  photons-per-bin  values  were  then  normalized  with  respect 
to  the  area  of  each  bin  lying  within  the  virtual  canopy’s  elliptical  footprint. 
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Figure  30  Simulated  Ground  Plane  Beam  Footprint  for  an  Illumination  Angle  of  80° 
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A  contour  plot  of  the  integrated  and  normalized  number  of  photons  arriving  within  each  bin  (i.e.,  the 
simulated  scattered  beam  footprint),  is  shown  in  Figure  31.  For  clarity  this  figure  has  been  scaled  to  a 
peak  value  of  unity  and  cropped  to  highlight  the  region  of  greatest  interest.  Note  that  all  photons  are 
initially  propagating  in  the  positive  range  (x)  direction  and  are  aimed  at  the  primary  target  located  at  the 
coordinate  system  origin  at  the  top  of  the  figure.  The  center  of  the  large  spot  on  the  ground  is  then 
positioned  beneath  the  location  where  photons  first  enter  the  canopy  (i.e.,  at  (x,y)  =  (-96.66,0)).  Recall 
that  the  majority  of  leaves,  facing  nominally  upwards,  will  be  distributed  near  the  top  of  the  canopy. 
Therefore  most  photons  are  scattered  downwards  just  after  entering  the  canopy,  which  in  turn  causes  the 
dominant  beam  footprint  to  be  found  beneath  the  entrance  location.  This  suggests  that  in  practice  it  might 
be  prudent  to  illuminate  the  canopy  directly  above  an  anticipated  target,  rather  than  at  an  oblique  angle. 
We  then  integrated  the  beam  footprint  data  across  both  the  range  and  cross-range  dimensions, 
individually,  and  calculated  a  best  fit  Gaussian  curve  to  describe  the  results.  Figure  32  contains  the  data 
(dots)  and  the  best  fit  Gaussian  curves  (solid)  for  both  the  range  (a)  and  cross-range  (b)  dimensions.  The 
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Figure  31  Simulated  Beam  Footprint  Cross  Sections,  and  Best  Fit  Gaussian  Distributions,  in  the  Range  (a) 

and  Cross-Range  (b)  Dimensions 
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Figure  32  Spatial  Distribution  of  Simulated  RMS  Pulse  Widths  on  the  Virtual  Canopy  Floor  (in  ns) 
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standard  deviations  of  the  best  fit  curves  (i.e.,  the  RMS  ground  spot  dimensions)  were  found  to  be  17.28m 
for  the  range  dimension  and  14.52m  for  the  cross-range  dimensions,  with  RMS  fitting  errors  of  9.1e-5  and 
2.0e-4,  respectively. 

Next,  we  created  a  contour  plot  of  the  RMS  temporal  pulse  width  (i.e.,  temporal  dispersion)  as  a  function 
of  spatial  location  as  shown  in  Figure  33,  where  the  gray  scale  legend  has  units  of  nanoseconds. 

Variations  in  path  length  lead  to  a  broad  spectrum  of  transit  times  through  the  canopy.  Therefore 
waveforms  measured  near  the  canopy  entrance  location  have  smaller  pulse  widths  than  do  those  measured 
further  away.  Notice,  for  example,  in  Figure  33,  that  the  dispersion  has  a  minimum  directly  beneath  the 
entrance  location  and  increases  as  the  distance  from  this  spot  increases. 

Additionally,  we  examined  the  temporal  data  for  photons  incident  upon  a  25  cm  diameter  virtual  detector 
pupil  located  in  the  receiver  plane.  The  1-D  arrival  time  probability  density  function  for  all  photons 
hitting  the  pupil  is  shown  in  Figure  34.  The  initial  shaip  rise  in  the  first  peak  is  due  to  the  influx  of 
photons  backscattered  from  high  within  the  canopy,  while  the  slowly  decaying  tail  results  from  multiply 
scattered  photons.  Photons  which  first  hit  the  secondary  target  under  the  canopy  entrance  location  also 
arrive  within  the  initial  peak,  but  are  inseparable  from  the  canopy  backscatter  noise.  Additionally,  there  is 
a  narrow,  shaip  peak  found  near  an  arrival  time  ofl. 7  microseconds  that  corresponds  to  those  photons 
which  propagate  to  and  from  the  primary  target  without  scattering  by  leaves.  These  “ballistic”  photons  are 
scattered  from  the  primary  target  back  toward  receiver  pupil  with  very  little  deviation  (i.e.,  a  very  nearly 
retro  reflection).  The  temporal  delay  between  the  initial  peak  and  the  unscattered  peak  exists  because 
photons  which  propagate  to  the  primary  target  and  back  travel  a  much  greater  round  trip  distance  than  do 
those  which  are  down- scattered  after  entering  the  canopy. 


Figure  33  Simulated  1-D  Temporal  PDF  Measured  in  Pupil  Plane  of  the  Virtual  Detector 


4.3  Experimental  Verification  Results 

A  sample  mean  waveform  captured  near  the  entrance  location  of  the  canopy  is  shown  in  Figure  36.  There 
are  two  large,  narrow  spikes  at  the  front  edge  of  the  pulse,  which  we  attribute  to  the  direct  reflection  from 
some  unidentified  scatterer,  followed  by  a  slowly  decaying  tail  as  was  seen  in  the  pulse,  Figure  36 
generated  in  our  simulation.  In  general,  these  captured  waveforms  have  more  structure  than  the  simulated 
waveforms,  presumably  due  to  the  clumping  and  clustering  of  leaves  and  branches.  Flowever,  the  general 
shape  of  the  measured  waveforms  is  consistent  with  that  of  the  simulation. 
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Figure  34  Image  of  the  Avalanche  Photo-Diode  Setup 


The  normalized  beam  footprint,  created  from  the  integrated  energy  contained  within  each  waveform,  is 
shown  in  Figure  37.  Note  that  all  photons  are  initially  propagating  in  the  positive  direction  along  the 
range  (x)  axis  and  are  aimed  at  the  primary  target  located  at  the  coordinate  system  origin.  As  predicted  by 
the  simulation,  there  is  a  spot  on  the  ground  centered  near  the  entrance  location  (i.e.,  at  (x,y)  =  (-76,0)) 
due  to  the  large  amount  of  down  scattering  just  after  photons  enter  the  canopy.  It  appears  that  there  was  a 
slight  misalignment  between  the  path  we  cleared  and  the  beam  trajectory,  however,  as  the  spot  on  the 
ground  is  seen  to  be  sloping  off  a  bit  to  the  right  side  in  the  plot.  We  attribute  this  to  our  cleared  path 
being  skewed  slightly  with  respect  to  the  actual  beam  trajectory.  Additionally,  there  appears  to  be  a  small 
discrepancy  in  our  illumination  angle,  as  we  expected  the  entrance  location  to  be  approximately  96.6m  in 
front  of  the  target.  At  the  primary  target  range  we  are  considering,  however,  this  corresponds  to  an 
illumination  angle  error  of  only  1.10°. 

We  then  integrated  the  beam  footprint  data  across  both  the  range  and  cross-range  dimensions, 
individually,  and  calculated  a  best  fit  Gaussian  curve  to  describe  the  results.  Figure  38  contains  the  data 
(dots)  and  the  best  fit  Gaussian  curves  (solid  lines)  for  both  the  range  (a)  and  cross-range  (b)  dimensions. 
The  standard  deviations  of  the  best  fit  curves  were  found  to  be  12.31  m  for  the  range  dimension  and  8.79 
m  for  the  cross-range  dimension,  with  RMS  errors  of  1.16e-4  and  5.01e-3,  respectively.  To  a  first  order 
this  data  matches  the  predicted  beam  footprint,  shown  in  Figure  31,  in  both  shape  and  size,  with 
discrepancies  likely  due  to  the  fact  that  our  simulated  data  represents  an  average  over  billions  of  canopy 
realizations,  whereas  our  experimental  data  arises  from  only  one  actual  canopy  realization. 
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Figure  35  Characteristic  Temporal  Waveform  Measured  by  an  Upwards  Looking  Avalanche  Photodiode 

Placed  on  the  Canopy  Floor 
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Figure  36  Measured  Ground  Plane  Beam  Footprint  for  an  Illumination  Angle  of  80° 
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Figure  37  Measured  Beam  Footprint  Cross-Section  Data,  and  Best  Fit  Gaussian  Curves,  in  the  Range  (a)  and 

Cross-Range  (b)  Dimensions 

Next,  the  temporal  dispersion,  or  RMS  pulse  widths  of  the  detected  signals,  is  shown  in  Figure  39.  where 
the  grayscale  legend  has  units  of  nanoseconds.  As  predicted  by  our  simulation,  the  measured  dispersion 
has  a  minimum  directly  beneath  the  entrance  location  and  increases  as  the  radial  distance  from  this  spot 
increases.  The  misalignment  between  our  cleared  the  path  and  the  beam  trajectory  can  also  be  seen  in  this 
figure.  Once  again,  discrepancies  between  our  simulated  and  experimental  dispersion  are  likely  due  to  the 
fact  that  our  simulated  data  represents  an  average  over  billions  of  canopy  realizations,  whereas  our 
experimental  data  arises  from  only  one  actual  canopy  realization. 

Finally,  from  the  tower  we  collected  a  series  of  5 1 1  range  gated  images  of  the  entire  illuminated  stand  of 
trees  using  a  gate  width  of  2ns  and  a  gate  delay  sequence  ranging  from  1.20ps  to  1.99ps.  For  each  frame, 
we  set  a  lower  threshold  to  reject  detector  noise  and  an  upper  threshold  to  eliminate  any  high  intensity 
first  surface  reflections.  Then,  by  integrating  the  pixel  values  within  each  thresholded  frame,  we  created  a 
one  dimensional  probability  density  function  for  the  photon  transit  time  through  the  canopy.  This  is 
shown  in  Figure  40.  The  shape  of  this  plot  is  very  similar  to  that  of  the  simulated  PDF,  shown  in  Figure 
33.  Both  are  characterized  by  an  initial  peak  due  to  the  high  density  of  leaves  near  the  top  of  the  canopy, 
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Figure  38  Spatial  Distribution  of  Actual  RMS  Pulse  Widths  Measured  on  the  Canopy  Floor 


Figure  39  Actual  1-D  Temporal  PDF  Measured  with  a  Range  Gated  Intensified  CCD  Camera  Located  in  the 

Tower 


each  with  a  width  of  approximately  100ns.  Following  the  initial  peak  is  a  period  of  very  low  return,  as  the 
few  photons  that  pass  through  the  top  of  the  canopy  experience  mostly  free  space  propagation  until  they 
reach  the  target.  Finally,  there  is  a  sharp  peak  around  1 ,75ps  corresponding  to  the  return  of  ballistic 
photons  which  strike  the  target  and  are  retro-reflected  back  towards  the  receiver. 


4.4  Monte  Carlo  Results  at  Various  Illumination  Angles 

Using  this  previous  work  as  a  basis  for  model  accuracy,  we  ran  our  Monte  Carlo  simulation  for 
illumination  angles  of  0°,  10°,  20°,  30  °,  60  °,  and  80 0  and  collected  the  spatial,  temporal  and  angular  data 
in  the  ground  plane  from  every  launched  photon.  We  then  divided  the  virtual  ground  plane  into  a  grid  of 
100x100  rectangular  bins  and  created  a  two-dimensional  histogram  of  the  number  of  simulated  photons 
arriving  in  each  bin  within  a  45°  half  field  of  view.  This  field  of  view  was  selected  to  match  that  of  the 
detector  used  in  our  experiment  described  later  in  this  section.  The  photons-per-bin  values  were  then 
normalized  with  respect  to  the  area  of  each  bin  lying  within  the  virtual  canopy’s  elliptical  footprint. 

A  contour  plot  of  the  integrated  and  normalized  number  of  photons  arriving  within  each  bin  (i.e.,  the 
simulated  scattered  beam  footprint)  is  shown  in  Figure  41  for  illumination  angles  of  0°  (a),  30°  (b),  60°  (c) 
and  80°  (d).  For  convenience  each  figure  has  been  scaled  to  a  peak  value  of  unity  and  cropped  to  highlight 
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the  region  of  greatest  interest.  Note  that  all  photons  are  initially  propagating  in  the  positive  range  (x) 
direction  and  are  aimed  at  the  line-of-sight  target  located  at  the  coordinate  system  origin  at  the  top  of  each 
figure.  The  center  of  the  large  spot  on  the  ground  is  then  positioned  approximately  beneath  the  location 
where  photons  first  enter  the  canopy  for  each  illumination  angle.  The  actual  photon  entrance  locations,  as 
determined  by  the  illumination  angle  geometry,  and  the  simulated  spot  centers  are  provided  in 


Table  7.  Recall,  in  our  model  the  majority  of  the  leaves,  facing  nominally  upwards,  will  be  distributed  near 
the  top  of  the  canopy.  This  in  turn  causes  the  dominant  beam  footprint  to  be  found  beneath  the  entrance 
location.  This  suggests  that  in  practice  it  might  be  prudent  to  illuminate  the  canopy  directly  above  an 
anticipated  target,  rather  than  at  an  oblique  angle. 

We  then  integrated  the  simulated  beam  footprint  data  across  both  the  range  and  cross-range  dimensions, 
individually,  and  calculated  best  fit  Gaussian  curves  to  describe  the  results.  Figure  42  shows  the  data 
(dots)  and  the  best  fit  Gaussian  curves  (solid)  in  both  the  range  (a)  and  cross-range  (b)  dimensions  for  the 
80°  data.  The  standard  deviations  of  the  best  fit  curves  (i.e.,  the  RMS  ground  spot  dimensions)  as  well  as 
the  RMS  fitting  errors  for  each  illumination  angle  are  then  provided  in 

Table  8. 
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Figure  40  Simulated  Beam  Footprint  on  the  Ground  for  0°  (a),  30 0  (b),  60 0  (c)  and  80 0  (d)  Canopy 
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Figure  41  Cross  Section  of  the  Simulated  Beam  Footprint  and  Best  Fit  Gaussian  Distribution  in  the  Range  (a) 

and  Cross-  Range  (b)  Dimensions 
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Table  7  Actual  Photon  Entrance  Locations  and  Simulated  Principle  Spot  Centers 


Ilium. 

Angle 

(deg) 

Range 
Entrance 
Location  (m) 

Range 

Spot 

Center  (m) 

0 

0 

0.0333 

10 

-3.174 

-3.644 

20 

-6.551 

-7.112 

30 

-10.392 

-10.80 

60 

-31.177 

-27.72 

80 

-102.08 

-96.66 

Table  8  Best  Fit  Gaussian  Beam  Widths  and  RMS  Fitting  Errors 


Ilium.  Angle 

(deg) 

Range 
Std  (m) 

Range 
Abs  Err 

Cross  Range 
Std  (m) 

Cross  Range 
Abs  Err 

0 

13.04 

1.608e-4 

12.96 

1.846e-4 

30 

14.08 

1.892e-4 

13.04 

1.690e-4 

60 

14.76 

6.366e-5 

13.80 

1.409e-4 

80 

15.52 

1.253e-4 

13.48 

2.037e-4 

Additionally,  we  created  one-dimensional  arrival  time  probability  density  functions  from  photons  incident 
upon  a  25cm  virtual  detector  pupil  located  in  the  receiver  plane,  as  shown  in  Figure  43  for  canopy 
illumination  angles  of  10°  (a)  and  30°  (b).  We  define  signal  photons  (gray  line)  to  be  those  which  strike 
the  target  positioned  beneath  the  entrance  location  and  noise  photons  (black  line)  as  those  which 
backscatter  from  the  canopy  without  target  interaction.  Flere  we  are  ignoring  the  photons  returning  from 
the  line-of-sight  target.  As  shown  previously  in  this  section,  illuminating  the  canopy  directly  above  an 
anticipated  target  projects  the  most  energy  onto  the  primary  target.  Had  we  known  this  prior  to  running 
our  extensive  simulation  we  may  not  have  included  a  target  in  line  with  the  initial  beam  trajectory.  Also, 
especially  at  higher  illumination  angles,  photons  that  hit  the  LOS  target  and  return  to  the  detector  plane 
are  almost  all  ballistic  photons;  that  is,  photons  that  propagate  through  the  gaps  between  the  leaves  to  and 
from  the  LOS  target.  These  ballistic  photons  represent  a  very  small  fraction  of  photons  returned  to  the 
detector  overall. 

In  practice  we  can  use  a  range  gated  ladar  system,  which  will  allow  us  to  open  and  close  the  shutter  on  the 
camera  at  very  specific  time  intervals.  The  time  between  creating  the  pulse  and  opening  the  shutter  on  the 
camera  is  referred  to  as  the  gate  delay,  and  the  length  of  time  the  shutter  remains  open  is  known  as  the 
gate  width.  By  using  a  gate  delay  approximately  equal  to  the  earliest  arrival  time  of  the  signal  photons  we 
may,  for  example,  ignore  the  initial  noise  peak,  thereby  increasing  the  signal-to-noise  ratio.  Also, 
adjusting  the  gate  width  will  allow  us  to  further  increase  the  SNR. 

Moreover,  multiple  leaf  interactions  not  only  lead  to  temporal  delays,  but  also  to  spatial  and  angular 
dispersion.  The  volume  scattered  noise  photons  are  typically  more  spatially  dispersed  in  the  canopy  than 
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are  signal  photons  and  therefore  generally  return  to  the  pupil  from  more  distant  locations  and  often  at 
larger  angles.  Thus,  we  are  able  to  further  mitigate  noise  by  only  accepting  photons  which  are  incident 
upon  our  detector  within  a  limited  field  of  view.  A  scatter  plot  of  the  arrival  time  versus  arrival  angle  for 
photons  incident  upon  the  detector  is  shown  in  Figure  44.  for  10°  (a)  and  30°  (b)  canopy  illumination 
angles.  Notice  that  many  of  the  noise  photons  (black  dots)  arrive  before  any  of  the  signal  photons  (gray 
dots)  and  those  that  arrive  at  the  same  time  generally  come  from  higher  up  in  the  canopy  and  therefore 
have  a  greater  angle  of  arrival. 

We  consequently  investigated  the  optimization  of  the  signal-to-noise  ratio  by  creating  two  three- 
dimensional  data  matrices  containing  the  number  of  signal  and  noise  photons  arriving  within  bins  of 
varying  gate  delay,  gate  width,  and  half  field  of  view.  The  gate  delays  spanned  the  range  from  0.2ps  to 
l.Ops  with  bin  spacings  of  O.Olps,  the  gate  widths  spanned  a  range  from  Ops  to  0.3ps  with  bin  spacing  of 
0.0  lps,  and  the  half  fields  of  view  spanned  a  range  from  0°  to  40°  in  2°  increments.  However,  rather  than 
calculating  the  pre-detection  optical  signal-to-noise  ratio,  found  by  dividing  the  number  of  signal  photons 
in  each  bin  by  the  number  of  noise  photons,  we  performed  a  more  in-depth  analysis  by  considering 
several  sources  of  post-detection  electrical  noise. 
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Figure  42  Photon  Arrival  Time  PDF’s  without  Range  Gating  or  Field  of  View  Filtering  for  10°  (a)  and  30“  (b) 
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Figure  43  Scatter  Plot  of  Photon  Arrival  Time  versus  Angle  of  Arrival  onto  Detector  for  10°  (a)  and  30°  (b) 


4.4.1.  SNR  Calculations 


In  addition  to  the  backscattered  photons,  there  are  several  sources  of  noise  which  are  inherent  to  the 
detection  process,  including  shot  noise,  thermal  noise,  and  dark  current  noise.  The  thermal  and  dark 
current  noises  are  function  of  the  detector,  and  are  constant  in  terms  of  the  number  of  detected  photons. 
Shot  noise,  however,  relies  heavily  upon  the  number  of  photons  collected  by  the  detector. 

The  post-detection  electrical  SNR  can  be  calculated  from  these  noise  sources  according  to  the  equation, 


SNR  = 


2 

dark  currents 


2 

signal 


2 

background 


2  ) 
thermal  / 


(37) 


This  equation  can  be  rewritten  in  terms  of  the  power  of  the  detected  photons  and  detector  responsivity,  R , 
according  to  the  following  relation, 


SNR  = 


_ ^(p,T)) _ 

2qB(lD  +fl(P,(0>)+R2{P,K0)  +  %^ 


(38) 


where  P s(t)  is  the  instantaneous  optical  signal  power,  P n(t)  is  the  instantaneous  background  optical  noise 
power,  and  P ,(t)  is  the  instantaneous  total  optical  power  (i.e.,  noise  +  signal)  [9,10].  Also  note  that  kB  is 
Boltzmann’s  constant,  T  is  the  temperature  (measured  in  Kelvin),  B  is  the  detector  bandwidth  (measured 
in  Hz),  Rl  is  the  load  resistance,  ID  is  the  dark  current  and  q  is  the  charge  of  an  electron. 
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For  a  given  gate  delay  and  half  field  of  view  the  instantaneous  optical  signal  power  can  be 
calculated  from  the  number  of  signal  photons,  Ns,  arriving  within  the  range  gate  according  to  the  equation, 


p,(<) 


At 


(39) 


The  mean-squared  optical  signal  power  can  be  found  by  integrating  this  expression  over  the  duration  of 
the  range  gate, 


+Tgw 

(P,2('))  =  ^r-  (p.T>*  (40) 

1  Sw  <o 

where  to  is  the  gate  delay  and  Tgw  is  the  gate  width.  We  can  approximate  this  value  by  breaking  the  range 
gate  into  M  temporal  bins  and  summing  the  number  of  signal  photons,  Nsj,  arriving  within  each  time  bin 
according  to  the  equation, 


1  M 


hvN. 


V  , 


A  t 


After  simplification,  the  mean-squared  optical  signal  power  can  be  expressed  as, 


(41) 


ihvY 

TgwAf 


(42) 


Finally,  by  applying  the  relation  At  =  Tgw  jM  and  the  scaling  parameter  a,  the  mean-squared  optical 
signal  power  can  be  expressed  as, 


M±(aNj  <43) 

1=1 

where  Tgw  is  the  gate  width,  h  is  Planck‘s  constant,  and  vis  the  optical  frequency.  The  total  number  of 
simulated  photons  launched  into  the  canopy  does  not  remotely  approach  the  number  of  photons  contained 
in  a  single  lOOmJ  pulse  from  a  1064nm  laser,  which  contains  about  5.4el7  photons  per  pulse.  In  order  to 
get  meaningful  results  we  needed  to  use  signal  and  noise  photon  values  that  would  be  on  the  order  of 
what  we  would  expect  for  a  typical  laser  pulse.  Therefore,  we  scaled  the  number  of  signal  and  noise 
photons  we  collected  in  our  simulation  by  a  factor  a,  the  ratio  of  the  typical  number  of  photons  per  pulse 
to  the  total  number  of  photons  we  have  sent  into  the  canopy.  The  scaling  parameters  we  used  in  this 
calculation  are  found  in  Table  12.  Note  that  a  varies  with  illumination  angle  as  we  only  simulated  enough 
photons  into  the  canopy  to  achieve  smooth  probability  density  functions,  and  because  a  larger  percentage 
of  launched  photons  returned  to  the  virtual  detector  at  lower  illumination  angles  than  at  larger  angles. 
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Similarly,  the  mean-squared  optical  background  power  can  be  found  from  the  number  of  background 
noise  photons,  N„j,  arriving  instantaneously  within  each  time  bin, 


(44) 


Finally,  the  mean  total  power  can  be  expressed  as  the  sum  of  the  signal  and  optical  background  power 
according  to 


PtU 


hv 


-a 


(N,+Nn) 


gw 


(45) 


where  Ns  and  Nn  are  the  total  number  of  signal  and  noise  photons,  respectively,  integrated  over  all  time 
bins. 


4.4.2.  Results  of  Field  of  View  and  Range  Gate  Filtering 

Consider,  now,  that  we  measure  the  return  of  photons  in  the  receiver  plane  with  a  PIN  photodiode  which 
has  the  following  parameters  at  a  wavelength  of  1064nm:  ID=  InA,  R  =  lOOmA/W  and  RL  =  1 000Q.  Let 
us  also  assume  a  bandwidth  of  B  =  2500Hz  and  a  temperature  of  T=  300K.  Then,  by  applying  these 
parameters  to  Eq.  (38)  and  using  the  number  of  signal  and  noise  photons  from  each  bin  in  our  3-D  data 
matrices,  we  can  construct  another  three-dimensional  data  matrix  containing  the  signal-to-noise  ratio  for 
each  HFOV/gate  width/gate  delay.  We  can  then  use  values  from  this  matrix  to  select  the  optimal  gate 
delay,  gate  width  and  half  field  of  view  for  achieving  the  best  SNR. 

A  contour  plot  of  the  SNR  as  a  function  of  both  half  field  of  view  and  gate  delay  for  a  constant  gate  width 
of  0.1  ps  is  shown  in  Figure  45  for  10°  (a)  and  30°  (b)  canopy  illumination  angles.  The  shape  of  these 
contour  plots  is  characteristic  of  the  SNR  as  a  function  of  half  field  of  view  and  gate  delay  for  most  gate 
widths.  In  general,  the  signal-to-noise  ratio  decreases  as  the  half  field  of  view  increases,  as  mostly  noise 
photons  are  returning  at  larger  angles.  Additionally,  the  SNR  rises,  peaks,  then  decays  as  the  gate  delay 
increases,  following  the  temporal  PDF  of  returning  signal  photons  seen  in  Figure  45.  More  importantly, 
there  is  a  small  pocket  near  a  gate  delay  of  0.4ps  where  the  signal-to-noise  ratio  is  relatively  large. 
Selecting  parameters  from  within  this  region  allows  us  to  achieve  a  large  boost  in  SNR.  Note  the  location 
of  this  high- SNR  region  shifts  to  larger  gate  delays  as  we  increase  the  illumination  angle  of  photons  into 
the  canopy  due  to  the  extra  optical  path  distance  added  by  moving  the  detector  further  from  the  target. 
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Table  9  Signal-to-Noise  Ratios,  Gate  Delays  and  Widths,  and  Half  Fields  of  View  for  Various  Canopy 

Illumination  Angles 


Illumination 
Angle  (deg) 

Scaling 

Parameter 

SNR  before 
filtering 
(dB) 

Gate 

Delay 

(ps) 

Gate 

Width 

(M-s) 

Half  field 
of 

View  (deg) 

SNR  after 
filtering  (dB) 

0 

5.65e9 

-11.39 

0.36 

0.02 

10 

8.88 

10 

3.15e9 

-11.79 

0.39 

0.06 

8 

6.61 

20 

1.51e9 

-11.74 

0.42 

0.10 

12 

2.52 

30 

4.02e8 

-12.98 

0.43 

0.04 

6 

-0.55 

60 

1.23e8 

-20.77 

0.54 

0.06 

8 

-0.90 

80 

3.04e7 

-29.57 

1.44 

0.10 

10 

-15.20 

We  located  the  optimum  signal-to-noise  ratio  from  well-populated  bins  within  this  pocket.  It  is  important 
to  note  that  optimum  does  not  necessarily  mean  maximum,  as  there  are  several  bins  that  are  populated  by 
only  a  few  photons.  The  SNR  values  for  these  bins  were  specific  to  the  exact  realization  of  our  simulation 
and  not  representative  of  the  actual  signal-to-noise  ratio.  We  therefore  disregarded  angular/temporal  bins 
with  fewer  than  50  returning  photons. 


We  found  that  using  a  0.39ps  gate  delay,  a  0.06ps  gate  width  and  an  8°  half  field  of  view  produced 
optimal  signal-to-noise  ratio  improvement  in  the  10°  illumination  angle  data,  while  a  0.43ps  gate  delay, 
0.04ps  gate  width,  and  a  6°  half  field  of  view  optimally  enhanced  the  SNR  of  the  30°  data.  We  calculated 
that  applying  these  filters  increased  the  signal-to-noise  ratio  of  the  10°  and  30°  data  from  -111  .39dB  to 
8.88dB  and  -12.98dB  to  -0.55dB,  respectively.  Figure  46  displays  the  one-dimensional  arrival  time 
probability  density  functions  previously  seen  in  Figure  43  after  both  range  gating  and  field  of  view 
filtering  have  been  applied.  Note  that  the  signal  is  now  larger  than  the  noise  for  the  10°  data  and  on  the 
same  order  as  the  noise  for  the  30°  data.  The  calculated  SNR  values,  optimum  range  gates  and  half  field 
of  view  filters  for  several  other  illumination  angles  are  found  in  Table  9  with  generally  the  same  order  of 
increase  in  SNR. 
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Figure  44  Contour  Plot  of  SNR  as  a  Function  of  Gate  Delay  and  Half  Field  of  View  after  Field  of  View  for  a 
Gate  Width  of  O.lps  for  10°  (a)  and  30°  (b)  Canopy  Illumination 
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(a) 


(b) 


Figure  45  Photon  Arrival  Time  PDF’s  after  both  Range  Gate  and  Field  of  View  Filtering  for  10°  (a)  and  30° 

(b)  Canopy  Illumination 


4.4.3.  Additional  SNR  Considerations 

While  we  included  the  total  average  photon  count  due  to  background  as  a  noise  term  in  the  denominator 
in  Eq.  (38),  there  are  several  instances,  for  example  thermal  imaging,  in  which  only  the  shot  noise  due  to 
the  background  is  included  in  the  denominator. 


SNR  = 


2  qB(lD  +  *<P,  (<)))+ 

V  )  . 


(46) 


To  illustrate  the  effect  of  only  keeping  the  background  shot  noise  contribution,  we  recalculate  Table  9  to 
show  the  SNR  for  the  case  that  the  total  average  photon  count  due  to  background  is  removed  from  the 
denominator.  This  is  shown  in  Table  10. 

Table  10  Signal-to-Noise  Ratios,  Gate  Delays  and  Widths,  and  Flalf  Fields  of  View  for  Various 

Canopy  Illumination  Angles 


Ilium. 

Angle 

(deg) 

Sealing 

Parameter 

SNR  before 
Filtering 
(dB) 

Gate 

Delay 

(ps) 

Gate 

Width 

(US) 

Half  field 

Of  View 
(deg) 

SNR  after 
filtering 
(dB) 

0 

5.65e9 

-11.39 

0.36 

0.02 

10 

15.95 

10 

3.15e9 

-11.79 

0.39 

0.06 

8 

13.23 

20 

1.51e9 

-11.74 

0.42 

0.10 

12 

7.61 

30 

4.02e8 

-12.98 

0.43 

0.04 

6 

3.47 

60 

1.23e8 

-20.77 

0.54 

0.06 

8 

-0.89 

80 

3.04e7 

-29.57 

1.44 

0.10 

10 

-9.46 
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5.  Conclusions  and  Recommendations 


We  have  presented  a  Monte  Carlo  model  used  to  simulate  the  propagation  and  scattering  of  light  through 
a  dense  tree  canopy.  The  model  is  characterized  by  the  probability  density  functions  governing  several 
parameters,  such  as  the  leaf  number  density,  the  angular  orientation  of  leaves,  and  the  bidirectional 
scattering  distribution  functions  of  individual  leaves.  We  then  measured  the  physical  dimensions  and  the 
leaf  area  density  of  a  nearby  grove  of  trees,  as  well  as  the  angle  and  distance  between  our  tower  and  the 
target  location  within  this  tree  grove.  We  then  applied  these  values  to  the  model. 

We  ran  our  simulation  for  80°  illumination  and  examined  the  expected  beam  footprint  and  pulse  width 
profile  on  the  ground,  as  well  as  the  temporal  returns  we  would  expect  to  measure  at  the  receiver.  Then,  in 
order  to  validate  the  model,  we  illuminated  our  tree  grove  with  a  780  nm  beam.  We  measured  waveforms 
on  the  canopy  floor  with  a  5x12  grid  of  upwards  facing  avalanche  photodiodes  and  collected  range  gated 
images  using  an  ICCD  camera  located  in  the  tower. 

As  predicted  by  the  simulation,  we  experimentally  verified  that  a  large  number  of  photons  are  down 
scattered  just  after  entering  the  canopy,  creating  a  large  spot  on  the  ground  beneath  the  entrance  location. 
We  then  found  that  a  Gaussian  distribution  fit  the  range  and  cross-range  cross-sections  of  the  simulated 
and  measured  beam  footprints  with  great  accuracy.  We  also  found  that  the  standard  deviations  of  the  best- 
fit  Gaussian  distributions  were  on  the  same  order  in  both  simulation  and  experiment.  Additionally,  we 
observed  that  the  pulse  widths  of  the  waveforms  measured  on  the  ground  were  similar  in  both  shape  and 
magnitude  to  those  predicted  by  our  simulation.  Finally,  we  found  that  the  measured  1-D  temporal  PDF 
of  photons  returning  to  the  receiver  closely  matched  the  PDF  predicted  by  the  simulation. 

We  believe  that  any  mismatch  between  our  simulation  and  experimental  results  can  be  attributed  to  two 
factors.  First,  in  the  experiment  there  was  a  slight  misalignment  between  the  trajectory  of  the  initial  beam 
axis  and  the  path  cleared  in  the  tree  grove,  which  resulted  in  a  small  angular  error  in  the  pointing  of  our 
beam.  Second,  and  most  importantly,  the  simulation  results  represent  an  average  over  billions  of 
realizations  of  the  canopy  while  the  experiment  considers  only  one  real  forest  realization.  Considering 
these  factors,  we  believe  our  experimental  results  correlate  very  well  with  our  simulation,  leading  us  to 
conclude  that  our  model  is  valid. 

We  then  investigated  the  ability  of  the  canopy  propagation  model  to  handles  seasonal  canopy  variations. 
We  first  quantified  the  effects  of  seasonal  changes  on  two  forest  parameters:  the  bidirectional  scattering 
distribution  function  of  individual  leaves  and  the  maximum  leaf  area  density  of  the  canopy.  Then  by 
applying  these  parameters  to  our  Monte  Carlo  canopy  propagation  model  we  predicted  the  effect  of 
seasonal  changes  on  the  shape  and  size  of  the  beam  footprint  as  measured  on  the  ground.  We  found  that 
as  the  forest  health  declines  and  the  leaves  fall  from  the  trees  the  size  of  the  spot  on  the  ground  increases. 
We  then  experimentally  verified  the  results  of  the  model  by  illuminating  a  real  forest  with  a  780nm  beam 
and  sampling  the  beam  footprint  with  an  upwards  facing  APD  beneath  the  canopy.  The  same  trend  was 
seen  in  the  experimentally  measured  data  as  was  predicted  by  the  simulation. 

Finally,  we  used  our  Monte  Carlo  canopy  propagation  model  to  make  predictions  about  detecting  objects 
through  the  canopy  at  several  different  illumination  angles.  We  have  used  our  Monte  Carlo  propagation 
model  to  simulate  the  scattering  of  photons  through  a  dense  canopy  using  several  different  illumination 
angles.  We  found  that  the  scattered  beam  footprint  on  the  canopy  floor  is  primarily  located  beneath  the 
entrance  location  of  photons  into  the  canopy  for  all  illumination  angles.  As  a  result,  illuminating  the 
canopy  above  an  anticipated  target,  rather  than  aiming  directly  at  it,  places  more  photons  on  the  target  and 
immensely  increases  the  signal-to-noise  ratio  of  detectable  returning  photons. 
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Additionally,  we  examined  the  effects  of  using  a  narrow  field  of  view,  range  gated  camera.  We  found  that 
multiply-scattered  photons  spend  more  time  in  the  canopy  than  do  singly-scattered  photons,  and  therefore 
generally  arrive  at  the  virtual  receiver  plane  at  a  later  time.  The  use  of  a  range  gated  camera  allows  us  to 
look  past  the  initial  return  of  backscattered  noise  photons,  increasing  the  signal- to-noise  ratio.  Multiple 
canopy  scattering  also  causes  spatial  and  angular  dispersion  of  photons  throughout  the  canopy.  Thus, 
photons  arriving  at  later  times  most  likely  also  arrive  at  larger  angles,  which  allows  many  noise  photons 
to  be  filtered  using  a  narrow  field  of  view  camera. 

We  found  that  by  illuminating  directly  above  the  desired  target  and  using  a  narrow  field  of  view,  range 
gated  camera  we  found  that  we  can  boost  the  signal-to-noise  ratio  on  the  order  of  15dB.  While  our  results 
are  specific  to  the  exact  geometry  and  canopy  architecture  used  in  this  simulation,  we  have  demonstrated 
in  principle  that  we  may  boost  the  signal-to-noise  ratio  of  detected  multiply  scattered  photons  using  a 
narrow  field  of  view,  range  gated  camera,  ultimately  providing  the  possibility  of  imaging  obscured  targets 
embedded  within  a  dense  forest  canopy  at  low  illumination  angles. 
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Appendix  A:  Calculation  of  Probability  of  a  Ballistic  Photon 


In  this  appendix  we  derive  the  expression  for  calculating  the  probability  of  encountering  an  unscattered 
photon  in  the  Monte  Carlo  simulation.  Before  proceeding  directly  into  the  derivation  for  calculating  the 
probability  of  an  unscattered  photon,  it  is  important  to  note  a  fundamental  statistical  rule  concerning 
conditional  probability.  Assume  that  A  and  B  denote  two  events  and  p(a\b)  denotes  the  conditional 

probability  of  A  occurring  given  that  B  has  already  occurred.  Baye’s  Theorem  then  gives  the  relationship 
between  the  two  stochastic  events  to  be, 


p(A\B)=p(B\A)^y 


(A.l) 


Recall,  in  the  simulation  we  have  broken  the  inhomogeneous  canopy  into  50  horizontal  slices  of  constant 
leaf  number  density.  We  did  this  so  that  the  probability  density  function  governing  the  random 
propagation  distance  within  each  region  could  be  modeled  as  a  negative  exponential,  expressed  as, 


pAd) 


1 

— exp 
D 


(A.2) 


where,  D  is  the  mean  free  path,  which  is  a  function  of  leaf  number  density. 


The  probability  that  a  photon  will  interact  with  a  leaf  within  a  region  \dx  ,d2  ]  of  constant  leaf 
number  density,  illustrated  in  Figure  61,  is  the  integral  of  the  PDF  over  that  region, 


d2 

P{dx  <  d  <  d2  }  =  |  pd  (d  )dd  . 

d , 


(A.3) 


Then  the  probability  that  the  photon  will  not  interact  with  a  leaf  in  the  region  [<:/  ,  d2  ]  is  equal  to  one 
minus  the  probability  that  there  will  be  an  interaction, 

p{dx  <d<d2}  = 

Because  of  the  segmented  canopy,  the  calculation  of  the  probability  of  an  unscattered  photon  must  be 
done  one  region  at  a  time.  So,  the  probability  of  not  having  an  interaction  in  the  ith  region  is  dependent 
upon  the  probabilities  of  there  being  an  interaction  in  any  of  the  previous  regions.  This  is  best  explained 
by  making  the  first  few  iterations  of  the  method. 

To  begin,  let  the  mean  free  path  of  the  ith  region  to  be  denoted  Di  and  the  boundaries  of  the  ith  region  to 
be  [dj_] ,  dt  ] .  Then,  the  probability  that  there  will  not  be  an  interaction  in  the  first  region  is  given  by  the 
expression, 


“2 

1  -  J pd(d)dd  . 


(A.4) 
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p,M) 


Figure  A-l  Probability  that  a  Photon  will  Interact  with  a  Leaf  within  the  Region  dl  to  d2  is  Given  by  the 
Integral  of  the  Negative  Exponential  PDF  over  that  egion 


PjO  <  d  <dl 


(A. 5) 


As  has  already  been  expressed,  the  probability  of  there  not  being  an  interaction  in  the  first  region  is  then 
given  by  one  minus  this  value, 


(A.6) 


The  probability  of  interaction  in  the  second  region  is  conditional,  given  that  the  photon  has  already 
propagated  out  of  the  first  region.  That  is,  let  A  denote  the  event  “photon  propagated  through  region  1 
without  interaction,”  and  let  B  denote  the  event  “interaction  occurs  in  region  2.”  Then  Baye’s  Theorem 
says  that  the  probability  that  the  photon  will  have  an  interaction  in  region  2,  given  that  it  has  already 
propagated  through  region  1  is  given  by, 


P(B\A)=  p[A\B)d^ 


(A.7) 


We  can  read  the  first  term  p[a\B  )  on  the  right  side  of  the  expression  to  be,  “the  probability  of  the  photon 

propagating  through  region  1,  given  that  an  interaction  occurs  in  region  2.”  If  the  interaction  occurs  in 
region  2,  then  the  photon  must  have  propagated  through  region  1  without  interaction.  Therefore,  this  value 
is  always  unity  and  Eq.  (B.7)  can  be  rewritten, 
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(A. 8) 


p(b\a)= 


Hb) 

P(A) 


To  evaluate  this  expression,  consider  Figure  A-2.  The  probability  that  an  interaction  occurs  in  region  2  is 
the  integral  of  the  PDF  over  the  boundaries  of  region  2, 


P{B }  =  P{dl  <  d  <d2}  = 


(A.9) 


This  region  is  illustrated  by  the  area  marked  with  negatively  sloped  lines.  The  probability  that  the  photon 
propagated  to  region  2  without  interaction  is  denoted  by  the  area  marked  with  the  positive  sloped  lines 
and  given  by  the  integral  of  the  PDF  from  the  beginning  of  region  2  to  the  extent  of  the  range, 
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(A.  10) 


Then,  the  probability  that  there  will  be  an  interaction  in  region  2,  given  that  the  photon  has  propagated 
through  region  1  is  found  by  substituting  equations  (A.9)  and  (A.  10)  into  expression  (A.8).  Subtracting 
this  from  unity  gives  the  probability  of  there  not  being  an  interaction  in  region  2,  given  that  there  was  no 
interaction  in  region  1 . 


(A.  11) 


p.M) 


Figure  A-2  The  Probability  Density  Function  that  a  Photon  will  have  an  Interaction 
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Finally,  one  can  extend  this  computational  stencil  to  the  ith  region. 


PR  !  <d<dt  dt_x  <  d  r=  1 


L  A 


f  1  d 

—  exp - dd 

j  n  n 


°°  1 

I  7Texp 


D, 


(A. 12) 


And  the  total  probability  that  there  will  not  be  any  interaction  throughout  the  entire  canopy  will  be  the 
product  of  all  the  probabilities  calculated  for  each  region, 


P\0<d  <h\=Y\  1 


f  1  d 

—  exp - dd 

D„  F  D, 


—  exp - dd 

D„  D. 


(A.13) 
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Appendix  B:  APD  Field  of  View  Characterization 


Prior  to  performing  the  waveform  measurement  experiments  we  characterized  the  field  of  view  of  the 
avalanche  photodiode  (APD).  Initially  we  had  planned  to  mount  a  narrow  band  pass  filter  on  the  APD  so 
that  we  could  perform  our  waveform  measurement  experiments  during  the  day,  however,  as  the  results 
within  this  appendix  show,  the  field  of  view  with  the  band  pass  filter  in  place  was  too  narrow. 

We  made  use  of  the  BRDF  measurement  apparatus  previously  used  to  measure  the  bidirectional  scattering 
distribution  functions  of  leaves,  which  is  depicted  in  Figure  B-l.  We  directed  a  linearly  polarized  (at  20° 
with  respect  to  vertical),  1064  nanometer  pulsed  laser  (>6pJ  pulse  energy)  along  the  axis  of  a  stationary 
optical  rail.  We  then  used  a  polarizer  and  two  wave  plates  to  attenuate  the  beam  to  avoid  damaging  the 
sensitive  APD. 


We  directed  the  reflection  from  the  second  wave  plate  towards  a  50%  reflective  Spectralon®  disc  and 
detected  the  pulses  with  a  InGaAs  PIN  diode.  We  used  the  measured  radiant  energy  reflected  from  this 
disc  as  both  a  trigger  to  signal  measurements  from  the  APD  and  also  as  a  pulse  energy  monitor  so  that 
data  could  later  be  normalized  with  respect  to  variations  in  pulse  energy. 

We  used  two  neutral  density  filters,  mounted  on  the  rail  after  the  beam  splitter  to  further  attenuate  the 
power  of  the  laser  in  order  to  avoid  damaging  the  APD.  The  ND  filters  were  tilted  slightly  in  order  to 
avoid  direct  reflections  into  the  beam  path.  Any  deflection  of  the  beam  due  to  the  first  filter  was 
counteracted  by  the  second  so  that  the  path  of  the  beam  remained  along  the  rail  axis. 


Si:  APD 
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Polarizer  + 
ND  Filters  Waveplates 


rff- - -<f4vr  Laser 
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Si:  PIN 

Figure  B-l  Schematic  Depiction  of  the  Experimental  Setup 


After  attenuation,  we  directed  the  laser  beam  onto  the  Hamamatsu  C5658  High  Speed  APD  Module 
which  was  mounted  on  a  motor  driven  and  computer  controlled  rotation  stage  which  also  tracked  the 
APD’s  rotational  position  relative  to  the  beam  path.  This  rotation  stage  is  capable  of  360°  of  rotation, 
allowing  the  APD  to  be  illuminated  at  any  angle. 

The  motors  of  both  the  detector  rail  and  leaf  mount  were  controlled  by  an  automated  LabView®  program 
which  rotated  the  detectors  and  leaf  to  predetermined  angles.  At  each  set  of  angles  the  program  collected 
512  samples  from  each  detector.  The  samples  were  then  averaged,  and  a  mean  value  was  stored  for  each 
detector.  The  program  then  moved  the  leaf  and/or  detectors  to  the  next  set  of  angles,  and  the  process  was 
repeated  until  all  the  desired  angle  permutations  are  exhausted. 
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Alignment  of  the  system  was  critical.  Slight  misalignment  of  the  active  area  of  the  APD  and  the  axis  of 
rotation  of  the  rotation  stage  would  result  in  very  large  errors.  We  aligned  the  system  by  mounting  a  cap 
with  crosshairs  onto  the  front  of  the  APD  such  that  the  crosshairs  were  centered  on  the  active  area.  A 
near-IR  camera  was  also  set  up  so  that  it  was  focused  onto  the  crosshairs/active  area  of  the  APD.  With  the 
laser  illuminating  the  crosshairs,  we  translated  the  rotation  stage  such  that  the  center  of  the  beam  fell  just 
above  the  “X”  on  the  crosshairs,  as  seen  in  Figure  B-2.  Because  the  crosshairs  were  mounted  out  in  front 
of  the  active  area  and  the  camera  was  mounted  above  the  laser  (looking  slightly  downward),  we  wanted 
the  beam  spot  to  be  slightly  higher  than  the  “X”  on  the  crosshairs.  This  placed  the  beam  directly  on  the 
center  of  the  active  area. 


Figure  B-2  The  Beam  is  Incident  on  the  Crosshairs  Just  Above  the  “X” 


We  placed  a  piece  of  scotch  tape  on  the  monitor  and  drew  dot  on  the  center  of  the  beam  spot,  which  can 
also  be  seen  in  Figure  B-2  above.  Then  we  turned  off  the  laser  and  positioned  the  APD  so  that  the  dot  on 
the  monitor  fell  exactly  on  the  center  of  the  active  area.  Note  that  only  the  two  transverse  directions  may 
be  aligned  in  this  fashion.  To  align  the  axial  component,  we  rotated  the  APD  to  some  angle,  say  30°,  then 
translated  the  APD  along  the  axis  until  the  dot  on  the  monitor  fell  exactly  at  the  center  of  the  active  area. 


Figure  B-3  Images  of  the  active  area  of  the  APD  (a)  at  0°  rotation  and  (b)  at  40°  rotation 
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Once  alignment  was  finished,  we  turned  off  the  lights  to  reduce  the  ambient  light  and  rotated  the  polarizer 
to  extinguish  as  much  of  the  beam  as  possible.  We  then  placed  ND  filters  with  optical  densities  of  400 
and  250  in  the  beam  path  to  further  diminish  the  signal.  Using  the  computer  and  the  “Interactive  Moving 
Leaves”  program  (LabView  code  written  by  John  Schmoll)  we  rotated  the  APD  to  the  desired  angles  and 
recorded  the  512  sample  average  at  each  angle. 

The  field  of  view  data  of  the  APD  with  the  band  pass  filter  in  place  is  found  in  Table  B-l  and  the  field  of 
view  data  of  the  APD  by  itself  is  found  in  Table  B-2.  The  field  of  view  of  the  detector  with  the  band  pass 
filter  in  is  only  on  the  order  of  10°,  whereas  that  of  the  APD  by  itself  is  approximately  50°.  This  is  in 
accordance  with  the  specifications  of  the  APD  module,  which  affirm  that  the  field  of  view  is 
approximately  48°. 


Table  B-l  Field  of  view  data  of  the  APD  with  the  band  pass  filter 


Angle 

(deg) 

APD 

Chi  Max  (mV) 

PIN  trigger 
Ch2  Max  (mV) 

0 

231.4 

31.98 

2 

225.3 

32.97 

4 

186.7 

32.67 

5 

115.9 

32.54 

6 

45.61 

32.50 

8 

0 

32.62 

-2 

225.6 

32.72 

-4 

212.9 

32.26 

-5 

201.7 

32.01 

-6 

144.8 

31.82 

-8 

22.20 

32.09 

-10 

0 

31.32 

Table  B-2  Field  of  view  data  of  the  APD  without  the  band  pass  filter 


Angle 

(deg) 

APD 

Chi  Max  (mV) 

PIN  trigger 
Ch2  Max  (mV) 

0 

438.9 

31.63 

5 

437.3 

31.7 

10 

415 

31.52 

15 

399.4 

31.47 

20 

370.4 

31.28 

25 

364.3 

32.41 

30 

322 

32.26 

35 

265.2 

32.54 

40 

224.4 

31.98 

45 

184.1 

32.12 

-5 

433.5 

31.67 

-10 

409.6 

31.27 

-15 

384.6 

31.59 

-20 

372.4 

31.57 

-25 

342.5 

31.35 

-30 

317.7 

31.99 

-35 

284.1 

30.69 

-40 

242.8 

31.07 
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The  data  in  the  absence  of  the  band  pass  filter  is  plotted  in  the  Figure  B-4.  The  drop  off  of  the  APD  signal 
as  a  function  of  angle  is  due  to  several  phenomena,  the  first  being  the  cosine  dependence  of  the  projected 
area.  As  the  APD  is  rotated  to  higher  angles,  less  of  the  active  area  is  projected  into  the  direction  of 
illumination.  The  angular  dependence  of  the  reflection  coefficient  at  dielectric  boundaries  also  adds  to  the 
dwindling  of  the  signal  as  the  incidence  angle  increases.  There  are  two  such  dielectric  interfaces:  (1) 
air/silicon  boundary  at  the  active  area  and  (2)  air/borosilicate  glass  boundary  at  the  window. 

Finally,  the  curve  seen  on  the  figure  is  the  theoretical  signal  expected  when  the  projected  area  and 
reflection  coefficients  are  taken  into  consideration.  The  theoretical  curve  is  normalized  to  the  maximum 
value  of  the  detected  signal,  which  happens  to  occur  at  an  angle  of  0°. 


Figure  B-4  Hamamatsu  C5658  APD  Module  Field  of  View  Characterization 
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Appendix  C:  Test  Canopy  GPS  Waypoints 


In  Section  3  we  discussed  measuring  the  dimensions  of  the  test  canopy  in  order  to  create  the  most 
reproducible  results  from  the  simulation.  The  data  and  details  of  the  calculation  are  found  in  this 
appendix. 

We  measured  the  dimensions  of  a  nearby  wooded  area  using  an  eTrex  Vista™  hand  held  GPS  unit.  We 
periodically  stored  waypoints  around  the  perimeter  of  the  tree  grove  and  converted  the  latitude/longitude 
coordinates  into  differential  distances  from  a  reference  location.  The  logged  waypoints  (solid  dots)  and 
corresponding  approximate  errors  (surrounding  circle)  for  each  point  are  illustrated  in  Figure  C-l.  Note 
that  the  reference  waypoint  is  marked  with  blue  in  the  following  figure.  The  data  for  the  locations  of  these 
waypoints  is  located  in  Table  C-l. 


Figure  C-l  The  Locations  of  the  Waypoints  are  Overlain  onto  a  Satellite  Image 


Conveniently,  the  tower  from  which  field  tests  are  conducted  is  directly  east  of  the  wooded  area. 
Therefore,  we  could  circumscribe  an  ellipse  around  the  woods  with  the  major  axis  equal  to  the 
longitudinal  distance  between  two  waypoints  on  the  horizontal  edge  of  the  canopy.  We  used  waypoints  22 
and  14  as  the  principle  locations  on  the  horizontal  axis.  The  two  are  separated  in  longitude  by  a  distance 
of  371  meters.  Similarly,  the  minor  axis  of  the  canopy  was  calculated  from  the  latitudinal  distance 
between  two  waypoints  on  the  vertical  edge  of  the  canopy.  We  used  waypoints  3  and  17  as  two  chief 
locations  on  the  vertical  axis.  The  two  are  separated  in  latitude  by  a  distance  of  229  meters.  The  ellipse 
circumscribing  the  tree  stand  is  shown  in  Figure  C-2. 
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Table  C-l  Longitudinal/latitudinal  locations  of  waypoints  and  differential  distance  from  the  reference 

waypoint** 


Waypoint 

Latitude 

Differential 
Distance  (m) 

Longitude 

Differential 
Distance  (m) 

46.674 

0 

5.1 

0 

2 

46.687 

24.01214376 

5.137 

68.34225533 

3 

46.714 

73.88351927 

5.189 

164.3908304 

4 

46.702 

51.71846349 

5.219 

219.8034698 

5 

46.692 

33.24758367 

5.257 

289.9928131 

6 

46.686 

22.16505578 

5.274 

321.3933088 

7 

46.692 

33.24758367 

5.287 

345.4054526 

8 

46.678 

7.388351927 

5.294 

358.3350685 

9 

46.706 

59.10681542 

5.285 

341.7112766 

10 

46.717 

79.42478322 

5.3 

369.4175964 

11 

46.694 

36.94175964 

5.322 

410.053532 

12 

46.676 

3.694175964 

5.325 

415.5947959 

13 

46.66 

-25.85923175 

5.327 

419.2889719 

14 

46.652 

-40.6359356 

5.311 

389.7355642 

15 

46.643 

-57.25972744 

5.289 

349.0996286 

16 

46.631 

-79.42478322 

5.287 

345.4054526 

17 

46.59 

-155.1553905 

5.268 

310.3107809 

18 

46.586 

-162.5437424 

5.206 

195.7913261 

19 

46.588 

-158.8495664 

5.17 

129.2961587 

20 

46.601 

-134.8374227 

5.152 

96.04857505 

21 

46.624 

-92.35439909 

5.109 

16.62379184 

22 

46.647 

-49.87137551 

5.11 

18.47087982 

Figure  C-2  The  Canopy  is  Taken  to  be  an  Elliptical  Cylinder  with  Minor  Axis  Length  of  115  Meters  and 

Major  Axis  Length  of  185  Meters 
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Appendix  D:  Cable  Delay  and  Dispersion 


In  Section  4  we  reported  on  the  experimental  verification  of  our  Monte  Carlo  canopy  propagation  model. 
In  these  sections  we  discussed  the  methodology  we  used  to  measure  the  waveforms  from  the  canopy 
floor.  Recall  that  we  sent  a  trigger  beam  towards  a  target  in  the  clearing  which  we  detected  with  a  PIN 
diode.  We  sent  the  signal  from  this  detector  through  a  500ft  cable  into  the  forest  into  the  oscilloscope, 
which  was  stationed  with  the  rest  of  our  data  collection  system.  We  then  sent  the  main  beam  into  the 
canopy  and  measured  the  waveforms  on  the  ground  with  the  APD,  which  was  also  connected  to  the 
oscilloscope. 

The  two  waveforms  collected  during  the  waveform  experiment  were  temporally  misaligned  due  to 
different  cable  lengths  connecting  the  detectors  to  the  scope.  The  objective  of  this  test  was  to  determine 
the  temporal  delays  so  that  the  waveform  data  may  be  synched. 

We  connected  the  output  of  a  function  generator  to  channel  1  of  an  oscilloscope  using  a  cable  of  arbitrary 
length  and  a  T-connector.  We  then  connected  one  end  of  the  500ft  trigger  cable  to  the  T-connector  and 
the  other  end  to  channel  2  of  the  scope.  We  then  generated  a  100kHz,  8ns  pulse  with  a  lOps  period  using 
the  function  generator  and  stored  the  pulses  after  propagation  through  the  cables.  The  amplitude  of  the 
pure  pulse  was  IV  peak-to-peak  and  the  edge  time  was  5ns. 

The  measured  pulses  using  the  500ft  trigger  cable  are  shown  in  Figure  D-l,  where  the  baseline  pulse  has  a 
maximum  at  3.6ns  and  the  trigger  pulse  has  a  maximum  at  795.2ns.  The  time  delay  between  the  pulses  is 
then  791.6ns.  Also,  the  baseline  pulse  width  is  9.4ns  and  the  trigger  pulse  width  is  22.2ns. 

The  measured  pulses  using  the  APD  cable  are  shown  in  Figure  D-2,  where  the  baseline  pulse  has  a 
maximum  at  3.6ns  and  the  ADP  pulse  has  a  maximum  at  59.2ns.  The  time  delay  caused  by  the  APD  cable 
is  then  55.6  ns.  Also,  the  APD  cable  broadened  the  pulse  width  from  9.4ns  to  9.6ns. 
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Appendix  E:  Selection  of  Random  Value 


Nearly  all  technical  computing  software  has  a  built-in  function  to  generate  uniformly  distributed  random 
variables  between  0  and  1 .  Often  times,  however,  we  desire  to  select  specific  values  of  a  random  variable 
described  by  a  non-uniform  probability  density  function.  The  procedure  for  doing  so  is  actually  quite 
simple. 

In  general,  consider  a  random  variable  Y  which  can  be  described  as  a  function  of  another  random  variable 
X  (i.e.,  Y  -  /(x)).  If  X  is  a  random  variable  uniform  on  [0,l) ,  and  if  we  know  the  probability  density 
function  pY{v )  for  Y ,  then  for  any  given  randomly  chosen  value  of  X  the  corresponding  random  value 
selection  for  Y  is  governed  by  the  relationship 


y  =  f{x)  =  FY\x). 

where  x  and  y  represent  specific  values  of  the  random  variables  X  and  Y ,  respectively,  and  where 
F  ( v) ,  known  as  the  Cumulative  Distribution  Function  (CDF)  of  Y ,  is  the  indefinite  integral  of  fY{y)  ■ 
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List  of  Acronyms,  Abbreviations  and  Symbols 


Acronym  Description 


APD 

BRDF 

BSDF 

BTDF 

CDF 

IR 

LAD 

ND 

NIR 

PDF 

RMS 

SAR 

SNR 

A 
T 
R 
0D 
0 , 
^INC 
0L 

A 

os 


0s 

9 

0 

0T 

0T 

0  G 


0G 

d 

D{z,6) 

A,(e) 

N(z) 


A 


0 


d  l 


Avalanche  Photo-Diode 

Bi-Directional  Reflection  Distribution  Function 
Bi-Directional  Scattering  Distribution  Function 
Bi-Directional  Transmission  Distribution  Function 
Cumulative  Distribution  Function 
Infrared 

Leaf  Area  Density 
Neutral  Density 
Near  Infrared 

Probability  Density  Function 
Root-Mean-Square 
Synthetic  Aperture  Radar 
Signal-to-Noise  Ratio 
Absorptance 
Transmittance 
Reflectance 
Detector  Angle 

Illumination  Angle 

Incidence  Angle 

Leaf  Zenith  Angle 

Leaf  Azimuth  Angle 

Leaf  Scattering  Zenith  Angle 

Leaf  Scattering  Azimuth  Angle 

Photon  Propagation  Zenith  Angle 
Photon  Propagation  Azimuth  Angle 
Target  Scattering  Zenith  Angle 
Target  Scattering  Azimuth  Angle 
Ground  Scattering  Zenith  Angle 
Ground  Scattering  Azimuth  Angle 

Photon  Propagation  Distance 
Mean  Free  Path 
Mean  Projected  Area 

Leaf  Number  Density 
Projected  Area 
Mean  Leaf  Area 
Mean  Leaf  Diameter 
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Acronym 

u(x  -  JC0) 

pAx) 

L(z ) 

R/n 

h 


R 

q 

B 

T 

Rl 


N,(t) 

Ps(t) 

N„(f) 

P»{>) 

h 


M 

a 

pR) 


Description 

Unit  Step  Function 
Probability  Density  Function 
Leaf  Area  Density 
Maximum  Leaf  Area  Density 
Canopy  Height 

Vertical  Canopy  Location  corresponding  to  Maximum  Leaf  Area  Density 

Detector  Responsivity 
Charge  of  Electron 
Detector  Bandwidth 
Temperature 
Load  Resistance 

Boltzmann’s  Constant 
Dark  Current 
Number  of  Signal  Photons 
Instantaneous  Optical  Signal  Power 
Number  of  Noise  Photons 
Instantaneous  Optical  Noise  Power 

Plank’s  Constant 
Frequency  of  Light 
Gate  Width 

Number  of  Temporal  Bins 
Photon  Scaling  Parameter 
Instantaneous  Optical  Detected  Power 
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